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Abstract 

This  thesis  explores  the  use  of  1,  and  mixed  optimization  methods  to  design 
flight  control  systems,  optimization  is  used  to  handle  tracking  issues  in  the  design  of 
digital  compensators.  Control  deflection  and  rate  limitations,  overshoot  and  undershoot 
limitations  and  steady-state  error  requirements  are  discussed.  Model-matching  techniques 
which  produce  acceptable  tracking  results  with  lower  order  controllers  are  also  examined. 
New  numerical  methods  for  continuous  Hj/Lj  and  discrete  optimization  are 
presented.  These  methods  are  used  to  design  an  aircraft  controller  in  continuous  and 
discrete  time  and  the  results  are  compared. 
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APPLICATIONS  OF  AND  MIXED  ll^/k  OPTIMIZATION 


L  Introduction 


1.1.  Overview 

Much  of  the  research  in  optimal  flight  control  design  in  recent  years  has  focused 
on  H2  andH^  optimization.  While  each  method  is  well  suited  for  specific  classes  of  inputs 
and  outputs  of  a  system,  neither  method  adequately  handles  the  complete  flight  control 
problem.  H2  optimization,  for  example,  attempts  to  minimize  the  energy  of  a  system’s 
output  to  a  white  Gaussian  noise  input.  The  resulting  H2  design  is  adept  at  handling 
noises,  but  can  have  poor  stability  margins.  Further,  additional  work  is  usually  required  to 
produce  an  H2  design  with  good  tracking.  optimization  attempts  to  minimize  the 
energy  of  a  system’s  output  to  an  unknown  but  bounded  energy  input.  The  resulting 
design  can  have  excellent  stability  margins  and  tracking  performance,  but  be  notably 
deficient  at  handling  noises.  Both  methods  are  extremely  poor  at  handling  “hard” 
magnitude  and  time  domain  constraints  on  the  system,  such  as  control  deflection 
limitations,  control  rate  limitations  and  overshoot  restrictions  in  the  system  response. 

The  optimization  method  minimizes  the  maximum  magnitude  of  a  system’s 
output  to  an  unknown  but  bounded  magnitude  input.  Since  this  optimization  method  is 
also  a  time  domain  approach,  it  can  handle  “hard”  magnitude  and  time  domain  constraints 
on  the  system.  Little  research  has  been  done  on  using  fj  optimization  to  improve  a 
system’s  stability  margins;  however,  it  has  been  shown  ([DDB94])  that  optimization  can 
produce  systems  with  good  tracking.  Unfortunately,  designs  can  be  deficient  at  handling 
noises. 

Although  optimization  still  requires  further  research,  it  is  clear  that  this 
optimization  method  alone  will  not  address  the  complete  flight  control  problem.  The  next 
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logical  step  to  a  complete  design  methodology  is  to  mix  some  of  the  optimization  methods 
discussed  above.  Many  researchers  have  worked  in  this  area,  but  few  have  incorporated 
optimization  into  a  mixed  design  approach.  A  notable  exception  is  Jacques  ([JR94]),  who 
has  developed  a  fixed-order  mixed  numerical  optimization  method  for  discrete 

systems.  This  approach,  however,  is  relatively  new  and,  thus,  questions  still  remain  on 
how  to  best  utilize  this  mixed  design  methodology.  A  mixed  design  approach  for 
continuous  systems  incorporating  the  continuous  version  of  the  problem,  known  as 
optimization,  has  not  been  developed. 

Given  that  some  mixed  optimization  method  is  desired,  the  next  question  is 
whether  this  method  should  be  utilized  in  continuous  or  discrete  time.  Even  though 
aircraft  are  inherently  continuous  systems,  most  aircraft  controllers  are  actually 
implemented  digitally.  One  way  of  producing  a  discrete  controller  in  this  case  is  to  use  a 
continuous  design  approach  and  discretize  the  resulting  controller.  Another  alternative  is 
to  discretize  the  continuous  plant  and  use  a  discrete  design  methodology  to  obtain  a 
discrete  controller  directly.  Perhaps  the  best  way  to  determine  which  approach  is  better  is 
to  produce  designs  using  each  approach  and  compare  the  results. 

This  thesis  explores  the  use  of  and  fixed-order,  mixed  optimization.  A 
fixed-order,  mixed  H^/Lj  design  approach  for  continuous  systems  is  also  developed  to 
compare  the  two  methods  of  producing  discrete  controllers  mentioned  above. 

1 .2 .  Review  of  Related  Work 

Dahleh  and  Diaz-Bobillo  ([DDB94])  have  done  the  most  comprehensive  work  on 
optimization.  They  pose  the  optimization  problem  as  a  linear  programming  problem  and 
solve  it  exactly  for  one  block  systems.  Three  methods  for  finding  approximate  solutions 
to  multi-block  problems  are  also  presented  in  [DDB94].  Dahleh  and  Diaz-Bobillo 
propose  a  few  methods  of  incorporating  control  deflection  limitations,  control  rate 
limitations  and  overshoot  restrictions  in  the  optimization  problem;  however,  many 
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implementation  details  are  omitted  and  few  comparisons  between  the  different  methods 
are  shown. 

Several  researchers  have  worked  on  fixed-order,  mixed-norm  optimization 
methods  for  continuous  systems,  but  Walker’s  approach  ([Wal94])  solves  the  most 
general  mixed-norm  problem.  [Wal94]  discusses  fixed-order,  mixed  optimization 
where  the  and  subproblems  can  be  posed  with  dissimilar  transfer  functions. 

Singular  and  multiple  constraints  can  also  be  incorporated  into  his  design  method. 
Walker  formulates  a  method  of  incorporating  Lj  constraints  into  his  algorithm  using  a 
forward  Euler  transformation  of  the  continuous  Lj  problem,  but  does  not  implement  the 
method. 

As  mentioned  in  the  previous  section,  Jacques  ([JR94])  has  developed  a  numerical, 
fixed-order,  mixed-norm  control  synthesis  method  for  discrete  systems.  Jacques’  work 
essentially  extends  Walker’s  method  to  discrete  time  systems,  and  incorporates 
optimization.  The  portion  of  Jacques’  algorithm,  however,  does  not  handle  systems 
with  both  fast  and  slow  dynamics  due  to  the  large  number  of  sample  periods  required  to 
estimate  the  pulse  responses  of  these  systems.  Since  most  aircraft  problems  have  fast  and 
slow  dynamics  (such  as  actuator  dynamics  and  phugoid  modes),  this  restriction  in  Jacques’ 
algorithm  severely  limits  the  number  and  type  of  aircraft  control  problems  that  can  be 
considered. 

1.3.  Research  Objectives 

The  purpose  of  this  research  is  first  to  thoroughly  investigate  and  compare 
different  magnitude  and  time  domain  constraints  that  can  be  added  to  the  optimization 
problem  to  produce  systems  with  good  tracking  characteristics.  Next,  fixed-order,  mixed 
H2/LJ  and  optimization  methods  will  be  developed  which  handle  realistic  flight 
control  problems  (i.e.  systems  with  both  fast  and  slow  dynamics).  Finally,  both  mixed 
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approaches  will  be  used  to  compare  two  different  methods  of  producing  discrete  aircraft 
controllers. 

1.4.  Outline 

This  thesis  consists  of  6  chapters  including  this  introduction.  Chapter  2  presents 
important  background  theory  that  is  used  throughout  this  thesis.  Topics  include  state- 
space  systems  and  transfer  functions,  transformations  from  continuous  to  discrete  time, 
the  Youla  parameterization,  linear  programming  and  duality  theory,  and'Lj  optimization 
methods,  and  Hj  optimization  (both  continuous  and  discrete). 

Chapter  3  discusses  using  optimization  to  solve  tracking  problems.  Unweighted 
and  weighted  sensitivity  minimization  problems  are  considered,  along  with  constraints  on 
the  control  deflections  and  rates,  steady-state  error,  and  maximum  overshoot  and 
undershoot  values.  Each  technique  is  demonstrated  on  a  realistic  flight  control  problem. 
Two  different  model  matching  methods  are  also  presented  which  yield  acceptable  tracking 
results  with  lower  order  controllers. 

Mixed  H^/L^  optimization  is  developed  in  Chapter  4.  New  numerical  methods  of 
evaluating  the  L;  norm  and  its  gradient  are  also  presented.  The  H2/LJ  algorithm  is  used  to 
design  a  series  of  continuous  aircraft  controllers  which  demonstrate  the  trade-offs  between 
pure  H2  and  Lj  designs. 

In  Chapter  5,  the  mixed  optimization  problem  is  presented.  New  numerical 
methods  of  evaluating  the  norm  and  its  gradient  are  developed  which  offer  dramatic 
improvements  in  terms  of  computational  efficiency  over  the  methods  presented  in  [JR94]. 
The  new  lyi^  algorithm  is  tested  on  a  discrete  version  of  the  aircraft  control  problem 
presented  in  Chapter  4,  and  several  mixed  designs  are  compared.  The  aircraft  model 
is  then  tested  with  a  discretized  controller  from  optimization,  and  the  results  are 
compared  to  a  similar  H^/l^  design.  Finally,  Chapter  6  presents  the  author’s  conclusions 
and  recommendations  for  further  research. 
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11.  Background  Theory 


2.1.  Chapter  Overview 

This  chapter  covers  some  of  the  basic  theory  used  to  conduct  the  research  in  this 
thesis.  State-space  systems  and  transfer  functions  are  discussed  in  Section  2.2.  Section 
2.3  covers  three  transformation  methods  from  continuous  to  discrete  time.  Section  2.4 
covers  the  Youla  parameterization.  Concepts  in  linear  programming  and  duality  which  lay 
the  foundation  for  discussing  fj  optimization  are  discussed  in  Section  2.5.  The  fj 
optimization  method  is  subsequently  developed  in  Section  2.6.  Section  2.7  covers  L, 
optimization.  The  chapter  ends  with  discussions  of  continuous  H2  optimization  in  Section 
2.8  and  discrete  optimization  in  Section  2.9. 

2.2.  State-Space  and  T ransfer  Functions 

All  finite  dimensional  linear  systems,  whether  continuous  or  discrete,  can  be 
written  in  state-space  form.  While  the  state-space  form  of  a  system  is  not  a  unique 
representation,  it  is  often  used  in  controller  synthesis  since  matrix  manipulations  are  easily 
handled  by  digital  computers. 

The  continuous  state-space  form  of  a  linear,  time-invariant  system  is  written  as 

x(t)  =  A,x(t)4-B,u(t) 
y(t)  =  C,x(t)4-D,u(t) 

and  the  discrete  state-space  model  is  written  as 

x(k-l-l)  =  A<iX(k)-l-BdU(k) 
y(k)  =  Qx(k)-f-D<jU(k). 

In  both  cases,  Ae  Be  91“^  ,  Ce  and  De  are  constant  matrices,  and  xe  91", 
ue  9tP  and  ye  91*^  are  the  state,  control  and  output  vectors,  respectively.  In  the  discrete 
state-space  model  x,  u  and  y  are  sequences  with  index  ke  {0,1,2,...}. 
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The  state-space  forms  above  are  considered  minimal  if  A  has  the  smallest  possible 
number  of  states.  A  state-space  system  is  minimal  if  and  only  if  (A,B)  is  controllable  and 
(A,C)  is  observable. 

The  Laplace  transform  is  used  to  transform  continuous  state-space  realizations  to 
input-output  transfer  functions.  The  z  and  A,  transforms  perform  the  same  function  for 
discrete  state-space  realizations.  Given  any  continuous  time  signal,  f,  starting  at  t=0  and 
sampled  at  T  intervals,  the  z  transform,  F(z) ,  is  defined  by 

F(z)  =  |;f(kT)z-\  (2.3) 

k=0 

where  z  =  e“^ .  Note  that  this  form  of  the  z  transform  can  only  be  used  on  causal  systems, 
i.e.  systems  with  current  outputs  that  are  not  influenced  by  future  inputs.  This  is  not  a 
significant  restriction  for  this  thesis,  however,  since  causality  will  be  assumed  for  all 
systems  considered.  The  A  transform,  F(A.) ,  of  a  continuous  signal  starting  at  t=0  is 
defined  as 

F(?l)  =  |;f(kT)A^  (2.4) 

k=0 

Using  the  Laplace  transform,  the  continuous  state-space  realization  can  be  put  in 
transfer  function  form, 

T^(s)  =  C,(sI-AJ-*B,+D,.  (2.5) 

The  z  transform  on  a  discrete  state-space  realization  produces 

T^(z)-C,(zI-A,)-'B,+D„  (2.6) 

and  the  A  transform  produces 

(\  V 

T^(A)  =  C,[^-I-A,J  B,+D,.  (2.7) 

Since  this  work  considers  both  continuous  and  discrete  systems,  great  care  is  taken 
to  differentiate  between  the  two  and  the  different  discrete  domains.  In  places  where 
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equations  or  results  apply  to  all  cases,  the  clarifying  notation  will  be  omitted.  In  addition, 
transfer  functions  in  this  thesis  will  occasionally  be  written  in  a  convenient  shorthand 
notation, 

=  C(^I-Ar‘B  +  D,  (2.8) 

where  ^=s,  z  or  1/X,. 

Control  system  design  is  standardly  done  in  terms  of  two  transfer  functions,  P  and 
K,  where  P  is  the  weighted  plant  and  K  is  the  compensator  to  be  found  (see  Figure  2.1). 


A 

B 

C 

D 

Note  from  Figure  2.1, 


m 

y 


=  p 


u  =  Ky 


where  m  is  a  vector  of  controlled  outputs  and  r  is  a  vector  of  exogenous  inputs, 
partitioned  as 


(2.9) 
If  P  is 


P  = 


(2.10) 


then  the  closed-loop  transfer  function,  ,  can  be  determined  from  the  lower  fractional 
transformation  (LFT) 
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(2.11) 


=  F,(P,K)  =  +P^„K(I-P^K)-^P^  . 

The  system  shown  in  Figure  2.1  is  well  posed  if  and  only  if  (I-P^^K)  is  invertible.  This 
condition  is  assumed  throughout  this  work. 

P  and  K  can  also  be  written  in  terms  of  the  state-space  matrices  as 


A 

B. 

B„ 

A 

B.' 

n 

D.. 

D. 

D.. 

K  = 

(2.12) 

c. 

LCk 

Since  the  system  is  assumed  to  be  well  posed,  both  (I-D^D^^)  and  (I-D^^DJ  are  invertible 
and  the  closed-loop  transfer  function  can  be  written  as 


(2.13) 


where 


A  = 


A-hB„(I-D,D^)-^D,C 


B.(I-D  a)-c 


B„(I-D,D^)-'C, 

-ll 


A,-hB,(I-D  D,)-'D  c 


(2.14) 


g  ^rB,  +  B„(I-D,D^)-^D,D^ 

^  [  B,(I-D^D,)-'D^ 

C™  =[C^-^D,,D,(I-D^DJ-‘C,  D„„(I-D,D^)-'C,] 

D.r=[D„.  +  D.uD.(I-DA)"'D^]- 


(2.15) 

(2.16) 
(2.17) 


Stabihty  for  both  continuous  and  discrete  systems  is  determined  by  the  eigenvalues 
of  the  A  matrix.  If  all  the  eigenvalues  of  a  continuous  system  lie  in  the  open  left-half  of 
the  complex  plane,  the  system  is  stable.  Discrete  systems  in  the  z-domain  must  have  all 
their  eigenvalues  inside  the  open  unit  disk  for  stability,  while  discrete  systems  in  the  X- 
domain  must  have  all  their  eigenvalues  outside  the  open  unit  disk.  A  distinct  relationship 
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exists  between  the  z  and  X  domains,  namely  z=l/X.  The  next  section  discusses 
relationships  between  the  continuous  and  discrete  systems. 

2.3.  Continuous  to  Discrete 

Most  often,  controller  designs  are  developed  for  continuous  systems  but  are 
actually  implemented  with  digital  controllers.  The  resulting  mixed  system  is  known  as  a 
sampled  data  or  a  hybrid  system.  One  effective  way  to  design  a  digital  controller  for  such 
a  system  is  to  discretize  the  system  and  complete  the  controller  design  in  the  discrete  time 
domain.  To  do  this  effectively,  however,  an  appropriate  method  must  be  used  to 
discretize  the  system. 

There  are  several  ways  to  transform  a  system  from  continuous  to  discrete.  Some 
transformations  maintain  the  integrity  of  the  stability  boundary,  others  do  not.  Some 
transformations  are  better  for  matching  frequency  responses  and  others  for  time  responses. 
Since  discrete  signals  only  contain  samples  of  their  continuous  counterparts,  all  the 
transformations  tend  to  distort  some  information.  The  three  transformations  that  are 
important  in  this  work  are  presented  below  with  their  corresponding  limitations. 

2.3.1.  Zero  Order  Hold  Equivalence 

In  actual  practice,  a  discrete  signal  does  not  directly  drive  a  continuous  plant.  The 
discrete  signal  is  passed  through  a  digital-to-analog  converter  (DAC)  which  produces  a 
continuous  output  from  the  discrete  input.  Most  DAC  devices  convert  a  binary  computer 
output  to  a  voltage  level  and  then  hold  the  level  for  T  seconds  until  the  next  computer 
output.  These  types  of  devices  are  known  as  zero  order  holds  (ZOH)  and  represent  the 
most  common  method  for  transforming  a  system  from  discrete  time  to  continuous  time. 

The  ZOH  method  is  also  known  as  the  step  invariant  method  because  it  matches 
the  step  response  of  a  discrete  system  to  a  continuous  one.  For  example,  the  step 
responses  of  the  continuous  system. 
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(2.18) 


Figure  2.2.  Step  responses  for  a  continuous  system  and  its  ZOH  equivalent 


The  main  limitation  of  ZOH  equivalent  systems  is  that  their  frequency  responses 
do  not  match  their  continuous  counterparts  well  unless  the  sample  rate  is  very  fast.  For 
this  reason,  the  ZOH  method  is  seldom  used  to  transform  continuous  filters  to  digital 
ones. 
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2.3.2.  Forward  Euler  Rule 


The  forward  Euler  rule  is  one  transformation  method  which  attempts  to  maintain 
the  frequency  response  of  the  continuous  transfer  function.  This  technique  essentially 
replaces  the  Laplace  variable,  s,  with  a  function  of  z.  The  mle  is  defined  as 

z  =  l  +  sT.  (2.20) 

Since  z  =  e^^ ,  it  is  easy  to  see  that  the  forward  Euler  rule  is  simply  the  first  two  terms  in 
the  Taylor  series  expansion  of  the  exponential. 

The  main  limitation  of  using  this  method  is  that  the  jco  axis  in  the  s  plane  does  not 
map  to  the  unit  circle  in  the  z-plane.  In  fact,  the  location  of  the  stability  boundary  is 
dependent  on  the  sample  rate  chosen  to  discretize  the  system.  Consider  the  system  in 
(2.18),  which  has  a  pole  at  -1.  The  forward  Euler  transformation  of  this  system  is 

H(z)  =  ^— .  (2.21) 

z-l  +  T 

which  has  a  pole  at  1-T.  Clearly,  if  the  sample  rate  chosen  is  greater  than  or  equal  to  2, 
the  resulting  discrete  system  will  be  unstable. 


2.3.3.  T ustin  or  Bilinear  T ransformation 

Like  the  forward  Euler  rule,  the  Tustin  or  bilinear  transformation  attempts  to 
maintain  the  frequency  response  of  the  continuous  transfer  function.  Unlike  the  forward 
rule,  however,  this  transformation  maintains  the  integrity  of  the  stability  boundary.  The 
transformation  is  defined  as 


z  =  ■ 


l+sT/2 


(2.22) 


l-sT/2 

To  see  that  the  stability  boundary  is  correctly  mapped  in  this  transformation,  let  s=jco  in 
(2.22)  and  solve  for  the  magnitude  and  phase  of  the  resulting  discrete  transfer  function 


l+jcoT/2 

l-joDT/2’ 


|zj  =  l,  0^  =2tancoT/2. 
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Notice  that  as  co  increases,  the  phase  increases,  but  the  magnitude  is  always  1.  Thus,  the 
jco  axis  in  the  s-domain  correctly  maps  to  the  unit  circle  in  the  z-domain. 

Since  the  Tustin  transformation  works  well  at  maintaining  frequency  and  stability 
information,  it  is  the  most  common  method  used  to  discretize  continuous  filters.  The  main 
disadvantages  of  using  this  method,  however,  are  that  it  does  not  maintain  the  step 
response  of  the  continuous  system  and  it  is  more  difficult  to  implement. 

A  comparison  of  frequency  responses  for  the  system  in  (2.18)  and  its  ZOH, 
forward  Euler  and  Tustin  transformations  with  T=0.5  sec  is  shown  in  Figure  2.3. 


Figure  2.3.  Frequency  response  comparisons 


Note  in  this  particular  example  that  the  sample  rate  is  fast  enough  for  the  ZOH  method  to 
match  the  frequency  response  of  the  continuous  system.  However,  the  sample  rate  is  too 
slow  for  the  forward  Euler  rule  to  provide  a  good  match. 
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2.4.  Youla  Parameterization 


Any  type  of  normed  optimization  on  the  closed-loop  transfer  function  described  in 
(2.13)  involves  a  search  over  all  possible  stabilizing  compensators.  This  search  can,  in 
general,  be  quite  difficult  since  the  closed  loop  transfer  function  is  a  nonlinear  function  of 
K  [see  (2.1 1)].  The  Youla  parameterization,  however,  defines  the  set  of  all  possible  real- 
rational  compensators  that  stabilize  the  closed-loop  transfer  function  in  terms  of  a  free 
parameter,  Q.  In  fact,  the  resulting  expression  for  K  is  affine  (linear  with  an  offset)  in  the 
parameter  Q.  Optimization  problems  can  thus  be  done  as  a  search  for  Q  instead  of  K, 
which  is  much  easier.  The  only  restriction  is  that  Q  must  be  a  stable,  real-rational  transfer 
function  (the  set  of  which  are  a  convex  set). 

To  understand  the  Youla  parameterization  further,  consider  the  system  in  Figure 
2.4,  where  K  =  F^(J,Q). 


Figure  2.4.  Youla  parameterization 


The  closed-loop  transfer  function  between  m  and  r  is  given  by 

T„,=F,(P,K)  =  F,(S,Q)  (2.23) 
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where 


S  =  F,(P,J). 


(2.24) 


If  S  is  partitioned  as 


S  = 


(2.25) 


(2.23)  becomes 

T..  =  S„+Sj,Q(I-S22Q)-’S,,.  (2.26) 

It  can  be  shown  ([Mac89])  that  the  Youla  parameterization  of  K  produces  822=0,  which 
simplifies  (2.26)  to 

T„,=S„+S„QS,„  (2.27) 

which,  for  ease  of  notation,  will  be  written  as 

T„,,  =  H-UQV 

from  this  point  on.  Note  that  this  expression  is  indeed  affine  in  terms  of  Q. 

It  can  also  be  shown  ([Mac89])  that  if  F  is  a  gain  matrix  which  stabilizes  A+B^^F  , 
and  L  is  a  gain  matrix  which  stabilizes  A+LCy,  then  H,  U  and  V  are  given  by 


■a+b„f  -B„F 

H  = 

0  A-fLCy 

b.+ld^ 

[  -D„„F 

(2.28) 


A-t-B„F 

K 

.C.+D,„F 

A-hLCy 

B.+LDy, 

.  Cy 

(2.29) 

(2.30) 


The  sizes  of  the  transfer  functions  U  and  V,  above,  are  often  used  to  classify 
different  types  of  problems.  If  U  and  V  have  the  same  number  of  rows  as  columns,  then  U 
and  V  are  square  and  all-pass,  or  inner.  Problems  of  this  type  are  known  as  one  block 
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problems.  If  U  has  fewer  columns  than  rows,  and/or  V  has  fewer  rows  than  columns,  the 
problem  is  classified  as  a  multi-block  problem. 

Note  that,  by  definition,  H ,  U  and  V  are  all  stable  transfer  functions.  Therefore, 
the  closed-loop  transfer  function  from  m  to  r  will  be  stable  if  and  only  if  the  parameter  Q  is 
stable. 

2 .5 .  Linear  Programming  and  Duality 

Since  linear  programming  and  duality  concepts  will  be  used  to  find  solutions  in 
the  next  section,  it  is  important  that  they  are  reviewed.  First,  consider  the  standard  linear 
programming  problem, 

min  C^X 

X 

subject  to  Ax>b  (2.31) 

Xi>0,  i  =  l,...,n 

where  x  is  a  column  vector  of  variables,  c  is  a  vector  of  coefficients,  and  A  and  b  are  a 
matrix  and  vector,  respectively,  which  form  constraints  on  the  variables.  Dimensionally,  x 
and  c€  SR",  be  SR'’  and  Ae  SR'””.  Any  optimization  problem  with  a  linear  objective  and  linear 
equality  and/or  inequality  constraints  can  be  transformed  into  this  form. 

For  each  standard  linear  minimization  problem  in  the  above  form,  known  as  the 
primal  problem,  there  corresponds  a  linear  maximization  problem  known  as  the  dual 
problem.  The  standard  dual  problem  is  written  as 

max  yb 

y 

subject  to  yA  <  c  (2.32) 

yj>0,  j=l,...,p 

where  ye  SR'’  is  a  row  vector  of  dual  variables.  Note  that  the  dual  problem  simply  reverses 
everything.  In  the  primal  problem,  c  was  the  cost  function  and  b  was  the  constraint.  In  the 
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dual  problem,  b  is  the  cost  function  and  c  is  the  constraint.  The  inequality  sign  also 
changes  between  the  two  problems,  and  the  unknown,  y,  becomes  a  row  vector  instead  of 
a  column  vector. 

The  primal  and  dual  problems  are  closely  related  through  several  key  results  of 
duality  theory.  First,  it  is  known  that  if  the  primal  or  the  dual  problem  has  an  optimal 
solution,  then  the  other  also  has  an  optimal  solution  and  their  values  are  the  same.  Second, 
if  X  and  y  are  any  feasible  vectors  in  the  minimum  and  maximum  problems,  then 

yb<c^x.  (2.33) 

This  concept  is  known  as  weak  duality.  Finally,  if  the  vectors  x  and  y  are  feasible  and 
c^x  =  yb ,  then  x  and  y  are  optimal. 

The  primal  and  dual  variables  are  also  related  through  the  complementary 
slackness  conditions.  If  x  and  y  are  feasible  solutions,  then  both  are  optimal  solutions  of 
their  respective  problems  if  and  only  if  for  all  i=l  ,...n, 

i)  Xi  >0=>(yA)i  =Ci, 

ii)  (yA)i<Ci=>Xi=0. 

These  conditions  are  also  known  as  the  alignment  conditions. 

2.6.  Optimization 

The  goal  of  optimization  is  to  minimize  the  maximum  magnitude  of  the 
controlled  output  of  a  system  given  a  bounded  magnitude  exogenous  input.  Vidyasagar 
([Vid86])  first  introduced  this  problem,  but  Dahleh  and  Pearson  ([DP87])  are  responsible 
for  its  more  general  solution.  Dahleh  and  Pearson’s  method  of  solution  involves  posing 
the  problem  as  a  linear  programming  problem.  The  goal  of  this  section  is  to  explain  then- 
method  of  solution. 
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To  simplify  the  explanation,  the  introductory  development  considers  the  case  of 
one-block  problems  only.  The  changes  necessary  to  find  solutions  to  multi-block 
problems  are  discussed  thereafter. 

The  system  in  Figure  2.1,  where  r(k)e9i”'  is  an  exogenous  input  sequence  of 
unknown  but  bounded  magnitude  and  m(k)e9i"“  is  the  output  sequence  to  be 
controlled,  represents  the  standard  problem.  If  <I)  is  the  closed-loop  transfer  function 
from  m  to  r,  then  the  objective  of  optimization  can  be  written  as 


min 

K  stabilizing 


min 

K  stabilizing 


max 


(2.34) 


Several  steps  must  be  taken  in  order  to  pose  this  as  a  linear  programming  problem. 
First,  the  nonlinear  absolute  value  function  in  the  norm  calculation  must  be  removed.  This 
is  accomplished  by  a  standard  change  of  variables  from  linear  programming.  Let 
<I>  =  -<h“ ,  where  and  are  sequences  of  n^  x  n^  matrices  with  positive  entries. 

The  norm  calculation  can  then  be  replaced  by 


+  (2.35) 

which  is  equal  to  the  norm  if,  for  every  (i,j,k),  either  (|)+  or  (j)'  is  zero.  But  either  (])+  or  (j)‘ 
must  be  zero  at  the  optimal  solution,  so  the  substitution  is  valid.  To  see  this  fact  more 
clearly,  consider  minimizing  the  sum  of  two  non-negative  numbers  that  must  remain  a 
fixed  distance  apart.  For  example,  consider  the  problem 

min(a-f-b) 

(2.36) 

subject  to  a  -  b  =  2,  a  >  0,  b  >  0. 

The  solution  to  this  problem  is  a=2,  b=0.  Note  that,  with  the  substitution  of  (j)+  and  (j)', 
(2.35)  is  linear  but  the  number  of  variables  has  doubled. 
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Before  searching  for  the  variables  and  <I>  which  minimize  the  one-norm  of  0, 
constraints  must  be  imposed  to  ensure  that  the  resulting  O  will  be  stable  and  realizable. 
These  two  problems  are  handled  by  the  Youla  parameterization.  Recall  that 

<I)  =  H-UQV,  (2.37) 

and  is  stable  if  and  only  if  Q  is  stable.  Also  recall  that  U  and  V  are  inner  for  one  block 
problems,  and  thus  invertible.  This  means  that  Q  can  be  solved  for  directly, 

Q  =  U''(H-d))V"\  (2.38) 

which  makes  it  easy  to  see  that  Q  will  be  stable  if  and  only  if  the  transfer  function  (H-O) 
cancels  the  unstable  zeros  of  U  and  V.  In  other  words,  if  the  unstable  zeros  of  U  and  V 
are  denoted  as  a.- ,  then  Q  will  be  stable  if  and  only  if 

d)(ai)  =  H(a.),  forl<i<N.  (2.39) 

Further,  if  d)  is  written  as  a  function  of  X, 

6(?i)  =  £d)(k)?i\  (2.40) 

k=0 

then  this  constraint  can  be  expressed  as 

'l  ai  af  •■•]rd>(0)l  rH(a,) 

1  a^  al  •••  4>(1)  _  H(a2) 

1  i  ;•  -  <i>(2)  "  : 

1  ajij  a^j  •••  :  _H(aj^) 

for  the  case  of  simple  zeros.  It  can  also  be  expressed  as  =  b^^^ ,  which  is  linear  in 

d). 

With  all  of  the  above  modifications,  the  optimization  problem  becomes 

min  2<i>^(k)  +  4)“(k) 

k=0 

subject  to  Af^^(<I)^(k)-<I)“(k))  =  bf,^^ 

0^(k)>0,<i>-(k)>0  (2.42) 
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which  is  a  linear  programming  problem  with  an  infinite  number  of  variables  and  a  finite 
number  of  constraints.  The  corresponding  dual  problem  has  a  finite  number  of  variables 
and  an  infinite  number  of  constraints.  However,  if  there  are  no  unstable  zeros  of  U  and  V 
on  the  unit  circle,  at  some  large  enough  k,  af  will  be  small  enough  that  only  a  finite 

number  of  these  constraints  will  be  active.  Thus,  the  dual  problem  is  finite  dimensional, 
and  an  exact  solution  can  be  found.  Further,  the  existence  of  a  solution  to  the  dual 
problem  guarantees  the  existence  of  the  same  solution  to  the  primal  problem.  In  fact,  the 
optimization  problem  can  be  solved  direcdy  in  the  primal  space  by  tmncating  the  above 

series  at  a  large  enough  value. 

There  are  a  few  modifications  that  must  be  made  to  the  above  formulation  for 
multi-block  problems.  First  of  all,  U  and  V  may  not  be  invertible.  However,  from  (2.29) 
and  (2.30)  it’s  clear  that  if  has  full  column  rank  and  has  full  row  rank,  then  a  left- 
inverse  of  U  and  a  right-inverse  of  V  will  exist,  which  is  all  that  is  required.  These  two 
restrictions  ensure  that  all  the  controls  are  penalized  and  that  no  measurements  are  perfect. 
Additionally,  for  multi-block  problems  it  is  the  left  unstable  zeros  of  U  and  the  right 
unstable  zeros  of  V  that  must  cancel  with  zeros  of  (H  -<1>) . 

Multi-block  problems  have  an  infinite  number  of  variables  as  well  as  an  infinite 
number  of  constraints,  and  thus  cannot  be  solved  exactly.  To  counter  this  problem, 
Dahleh  and  Diaz-Bobillo  ([DDB94])  proposed  three  ways  to  find  approximate  solutions. 
The  first  method,  known  as  the  Finitely  Many  Variables  (FMV)  approach,  constrains  the 
polynomial  solution  O  to  a  fixed  length.  The  resulting  compensator  provides  a  sub- 
optimal  but  feasible  solution  to  the  problem.  The  second  method,  known  as  the  Finitely 
Many  Equations  (FME)  approach,  truncates  the  number  of  dual  variables,  which  is  the 
same  as  solving  the  primal  problem  with  a  finite  number  of  constraints.  The  solution  to 
this  problem  is  super-optimal  and  infeasible.  The  final  and  most  viable  method  is  known  as 
the  Delay  Augmentation  (DA)  approximation.  This  method  is  generally  considered  the 
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best  method  to  use  for  multi-block  problems  since  it  carries  more  information  about  the 
optimal  solution  than  the  other  two  approaches. 

The  DA  approach  embeds  the  multi-block  problem  into  a  larger  one-block  problem 
by  augmenting  pure  delays  to  U  and  V.  The  resulting  one-block  problem,  which  contains 
extra  degrees  of  freedom  in  Q,  can  then  be  solved  exactly.  While  the  solution  to  this 
problem  is  super-optimal  and  infeasible,  it  serves  as  an  upper  bound  to  the  true  optimal. 
To  get  a  feasible  solution,  the  extra  degrees  of  freedom  are  simply  stripped  out  of  Q.  The 
resulting  solution  is  sub-optimal  but  provides  a  lower  bound  to  the  optimal  solution. 
Thus,  this  method  produces  both  a  feasible  solution  and  bounds  on  the  optimal  solution. 

The  optimization  method  is  discussed  further  in  Chapter  3.  In  parts  of  Chapter 
3,  the  standard  linear  programming  problem  in  (2.42)  is  augmented  with  constraints  on 

the  maximum  magnitude  of  the  controlled  output  to  an  exogenous  step  input.  These 
additional  constraints  are  posed  using  the  vector  norm, 

||m|L  =  sup|m(k)| .  (2.43) 

k 


2.7.  Lj  Optimization 

The  goal  of  Lj  optimization  is  the  same  as  fj  optimization,  except  that  it  is  applied 
to  continuous  systems.  While  optimization  attempts  to  minimize  the  absolute  sum  of  the 
pulse  response  of  a  discrete  system,  L,  optimization  attempts  to  minimize  the  absolute 
integral  of  the  impulse  response  of  a  continuous  system. 

The  impulse  response  matrix  for  a  continuous  causal  system. 


[A, 

B.1 

D.J 

T^(s)  = 


with  p  inputs  and  q  outputs  is  given  by 

H(t)  =  C„e"^^B,+D„,5(t), 


(2.44) 


(2.45) 
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where  6(t)  is  an  infinite  pulse  of  zero  time  duration  applied  at  t=0.  If  hjj  are  the  elements 
of  H(t),  then  the  objective  of  Lj  optimization  can  be  written  as 


min 

K  Stabilizing 


min 

K  stabilizing 


(2.46) 


Unlike  optimization,  there  is  no  known  way  of  solving  the  L,  problem  exactly. 
In  fact,  Dahleh  and  Diaz-Bobillo  ([DDB94])  have  shown  that  if  a  solution  could  be  found, 
the  optimal  compensator  would  have  infinite  order.  Theory  has  been  developed  to  find 
approximate  solutions  to  the  optimization  problem.  Much  of  this  theory  involves 
discretizing  the  continuous  system  with  a  forward  Euler  transformation  and  performing 
optimization  on  the  resulting  discrete  system.  It  can  be  shown  ([Wal94])  that  the  norm 
of  this  discrete  system  is  an  upper  bound  to  the  Lj  norm  of  the  original  continuous  system, 
and  that  the  norm  monotonically  approaches  the  L,  norm  as  the  sample  period  is 
decreased  to  zero. 


Unfortunately,  this  approach  is  not  very  practical  for  systems  with  slow  modes. 
Since  these  systems  generally  have  very  slowly  decaying  pulse  responses,  a  larger  sample 
period  is  required  to  capture  the  pulse  response  with  a  reasonable  number  of  samples. 
This  larger  sample  period,  however,  leads  to  poor  approximation  of  the  Lj  norm.  Further, 
it  may  not  be  possible  to  pick  a  large  sample  period  and  have  the  discrete  system  remain 
stable.  Recall  from  Section  2.3.2  that  the  forward  Euler  transformation  does  not  correctly 
map  the  stability  boundary  from  continuous  to  discrete  time,  and  larger  sample  periods 
increase  the  likelihood  that  the  discrete  system  will  be  unstable. 

Two  numerical  methods  of  computing  the  Lj  norm  of  a  continuous  system  are 
presented  in  Chapter  4.  Both  methods  avoid  the  above  complications  by  approaching  the 
problem  more  directly.  The  primary  focus  in  Chapter  4  is  fixed-order,  mixed  H2/L1 
optimization  for  continuous  systems.  Continuous  H2  optimization  is  discussed  in  the  next 
section. 


2-17 


2.8.  Continuous  Optimization 

The  goal  of  H2  optimization  is  to  find  the  internally  stabilizing  controller  which 
minimizes  the  energy  of  the  system  output  to  a  white  Gaussian  noise  input.  The 
continuous  time  development  is  shown  here,  and  the  discrete  time  version  is  shown  in  the 
next  section. 

Consider  the  standard  Hj  problem  in  Figure  2.5,  where  w  is  a  zero-mean  white 
Gaussian  noise  input  with  unit  intensity  and  z  is  the  controlled  output. 

Z 

y 

Figure  2.5.  H2  problem 


Given  this  setup,  the  Hj  optimization  objective  can  be  written  as 

ll^zwIL- 

K  stabilizing  ^ 


If  T„„  is  written  as 


then  the  two-norm  of  T^^  can  be  expressed  as 

||X.(s)g  =  tr[QClCJ, 

where  Q  is  the  positive  semidefinite  solution  to  the  Lyapunov  equation 

A,Q-fQA^-hBX=0. 


(2.47) 

(2.48) 

(2.49) 

(2.50) 
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P  can  also  be  written  in  terms  of  state-space  matrices  as. 


'a 

p= 

C. 

D.W 

Dzu 

>. 

u 

_ _ 1 

Dyw 

(2.51) 


The  following  assumptions  are  now  made  on  the  state-space  elements  of  P(s): 

i)  Dzw=0 

ii)  D^=0 

iii)  (A,  B  J  is  stabilizable  and  (Cy,  A)  is  detectable 

iv)  D  =  I  and  D  =  I 


A-jcoI 
C,  D,, 


A-j0)I  B, 
C  D 


has  full  column  rank  for  all  co 


has  full  row  rank  for  all  co 


Assumption  i)  is  required  in  order  to  ensure  that  will  have  a  finite  two-norm. 
Assumption  ii)  is  not  required,  but  it  simplifies  the  problem.  Assumption  iii)  is  required 
for  a  stabilizing  compensator  to  exist.  Assumption  iv)  ensures  that  and  Dy^D^ 

have  full  rank  which  keeps  the  problem  from  being  singular.  This  assumption  can  be 
relaxed  to  just  full  rank  requirements  on  and  Dy^Dj^^  through  scaling.  Finally, 

assumptions  v)  and  vi)  are  required  for  the  existence  of  stabihzing  solutions  to  the 
algebraic  Riccati  equations  used  in  the  H2  solution. 

The  unique  controller  which  minimizes  (2.47)  is  given  by 


K2  0.(s)  = 


- 1 

> 

K, 

1 - 

1 

0 

0 

(2.52) 


where 


Aj=A-K,Cy-B„K, 


(2.53) 
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k,  =  b:x+d^c. 


(2.54) 

(2.55) 


k,=yc;4-b„dJ;„. 

X  and  Y  are  the  real,  unique,  symmetric,  positive  semidefinite  solutions  to  the  algebraic 
Riccati  equations 

(A-B„D^CJ"X  +  X(A-B„DLCJ-XBXX  +  C:C,  =0  (2.56) 

(A-B^Dj,C,)"Y  +  Y(A-B,D;,q)-YC,CjY  +  B:B,  =0,  (2.57) 

where 

C,  =  (I-D^Dl)C,  (2.58) 

and 

B„  =  BJI-Dj:„D,J.  (2.59) 

The  minimum  two-norm  which  corresponds  to  (s)  is  denoted  as  a . 

2.9.  Discrete  H2  Optimization 

While  the  goal  of  optimization  does  not  change  in  the  discrete  time  case,  there 
are  some  distinct  differences  between  the  discrete  and  continuous  time  solutions.  The 
purpose  of  this  section  is  to  highlight  these  differences. 

T^(z)  can  also  be  expressed  in  terms  of  the  state-space  matrices  in  (2.48).  The 
two-norm  of  T^^  in  discrete  time  is  given  by 

||T„(zf  =  tr[D„Dl,  +C.QCI],  (2.60) 

where  Q  is  positive  semidefinite  solution  to  the  discrete  Lyapunov  equation 

A3QA^-hBX=Q.  (2.61) 

Given  the  state-space  description  of  P(z)  in  (2.51),  let 

St  D.]  =  t  t 

L^zuJ  L^C  *<-cJ 


2-20 


and 


1 

PQ 

i _ 

[b'^ 

]  = 

r 

a 

Sf‘ 

L  w 

Ls; 

Rf. 

(2.63) 


The  following  assumptions  are  now  made  on  the  state-space  elements  of  the  plant,  P(z), 
using  the  above  definitions: 

i)  (A,  BJ  is  stabilizable  and  (C^,,  A)  is  detectable 

ii)  Rf,  R^>0 

iii)  D^C,  =  0and  B„Dj„=0 


Assumption  i)  is  required  in  order  for  a  stabihzing  compensator  to  exist.  Assumption  ii) 
ensures  the  existence  of  stabilizing  solutions  to  the  discrete  algebraic  Riccati  equations. 
Assumption  iii)  is  a  standard  assumption  which  is  not  required  but  simplifies  the 
derivation. 

The  unique  discrete  H2  optimal  compensator  is  given  by 


where 


Kaop.(z)  = 


A-AK,Cy-B„K,K,C^ 

ak,-b„k,k; 

-K,-FK,KfCy 

-k,k; 

K^=iR^  +  BlXBJ-\Sl  +  BlXA) 

K,  =  (YCJ;  -fSj)(C,YCj H-R,  +  C,S^  +S,C]r^ . 


(2.64) 


(2.65) 

(2.66) 


X  and  Y  are  the  real,  unique,  symmetric,  positive  semidefinite  solutions  to  the  discrete 


algebraic  Riccati  equations 

A^XA  Q,  -  ( A^XB„  +  S, )(R,  -h  B^XB„  )-^  (B^XA  -f-  S^)  =  X  (2.67) 

AYA"^  +  Qf  -  (AYCj  H-  Sj )(Rf  +  YCj)“'  (C^ YA^  +  )  =  Y  .  (2.68) 


Notice  from  (2.52)  and  (2.64)  that  Dj^  equals  zero  in  the  continuous  time  case  but 
it  is  not  equal  to  zero  in  the  discrete  time  case.  This  fact  is  specifically  pointed  out 
because  many  references  in  the  control  literature  incorrectly  state  that  Dj^  equals  zero  in 
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the  discrete  solution.  In  fact,  assuming  that  equals  zero  generally  results  in  a  sub- 
optimal  solution. 

2.10.  Chapter  Summary 

This  chapter  developed  some  essential  background  theory  which  will  be  required  in 
Chapters  3-5.  The  theory  on  L;  and  continuous  H2  optimization  will  be  important  in 
Chapter  4,  where  mixed  H2/LJ  optimization  for  continuous  systems  is  discussed.  Discrete 
H2  optimization  theory  is  used  in  Chapter  5  to  develop  a  mixed  optimization  method 
for  discrete  systems.  The  next  chapter  discusses  using  optimization  to  solve  tracking 
problems. 
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III.  Using  lx  Optimization 


3.1.  Chapter  Overview 

This  chapter  focuses  on  using  1^  optimization  to  solve  tracking  problems.  Section 
3.2  describes  some  changes  that  had  to  be  made  to  the  optimization  software  prior  to 
conducting  most  of  the  work  in  this  chapter.  Section  3.3  discusses  optimization  on  a 
conceptual  level,  and  brings  up  an  important  relationship  between  the  unit  pulse  and  step 
responses.  Unweighted  and  weighted  sensitivity  minimization  are  discussed  in  Section 
3.4.  Sections  3. 5-3.7  cover  different  constraints  that  can  be  added  to  the  optimization 
problem  to  handle  several  tracking  issues.  The  last  section  in  this  chapter  discusses  the 
use  of  ?!  optimization  for  model  matching  problems. 

Throughout  the  chapter,  a  Single  Input  Single  Output  (SISO)  longitudinal  model 
of  the  AFTI  F-16,  shown  in  detail  in  Appendix  A,  is  used  to  illustrate  various  tracking 
design  issues.  The  tracking  problem  described  in  this  chapter  is  defined  as  the  ability  to 
accurately  command  a  Ig  normal  acceleration  of  the  aircraft.  The  stabilator  is  the  only 
control  surface  considered  in  the  model,  and  it  is  limited  to  ±  25  degrees  deflection  angle 
and  ±60  degrees/sec  deflection  rate.  In  Section  3.8  a  model  matching  design  is 
completed  for  this  system  and  also  for  a  Multiple  Input  Multiple  Output  (MIMO)  system 
involving  a  missile.  The  objective  of  the  missile  problem  is  stated  in  Section  3.8  and  a 
detailed  description  of  this  system  is  given  in  Appendix  B. 

All  simulations  in  this  chapter  are  done  with  sampled-data  systems  (i.e.  discrete 
controllers  with  continuous  system  models).  Step  inputs,  applied  one  second  after 
simulations  are  started,  are  used  to  evaluate  tracking  performance.  Plots  of  the  control 
rate  are  based  on  finite  difference  calculations  over  the  discretization  period  of  the  system. 
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3.2.  fj  Optimization  Software 

MATLAB'^'^  software  used  to  perform  (j  optimization  was  first  written  by  Diaz- 
Bobillo  in  1991-92.  This  research  code  was  used  to  run  the  examples  in  the  book,  Control 
of  Uncertain  Systems:  A  Linear  Programming  Approach,  written  by  Dahleh  and  Diaz- 
Bobillo  ([DDB94]).  A  copy  of  the  software  and  a  preprint  of  the  book  were  used  to 
conduct  the  research  in  this  thesis. 

After  running  several  examples  with  Diaz-Bobillo’s  software,  it  became  apparent 
that  certain  “bugs”  existed.  The  first  error  occurred  with  certain  simple  one-block 
systems.  While  the  main  routine  always  displays  a  lower  and  upper  bound  to  the  one- 
norm  for  all  systems,  the  two  bounds  should  be  the  same  for  one-block  systems  since  they 
can  be  solved  exactly.  Due  to  numerical  differences  in  the  ways  the  upper  and  lower 
bounds  are  calculated,  the  two  bounds  are  never  exactly  the  same.  However,  for  several 
one-block  systems  there  were  wide  disparities  between  the  two  values.  The  problem  was 
traced  to  a  routine  which  rounded  off  the  unstable  zero  frequencies  of  U  and  V.  The 
amount  of  round-off  affected  the  interpolation  conditions  given  in  (2.41)  and  thus  the 
norm  calculations.  This  problem  was  corrected  by  eliminating  the  round-off  algorithm 
altogether. 

The  second  problem  encountered  with  Diaz-Bobillo’s  code  was  that  it  did  not 
directly  ensure  the  required  existence  of  a  left  inverse  of  U  and  a  right  inverse  of  V.  It 
appears  that  Diaz-Bobillo  attempted  to  address  this  issue  by  changing  the  discrete  system 
to  continuous  time  through  a  bilinear  transformation.  In  many  cases,  this  action  leads  to  a 
system  where  D^^,  has  full  column  rank  and  D^  has  full  row  rank,  which  ensures  that  U  has 
a  left  inverse  and  V  has  a  right  inverse.  Unfortunately,  these  conditions  must  be  met  for  a 
discrete  time  system,  not  a  continuous  one,  and  a  bilinear  transformation  back  to  discrete 
time  recreates  the  rank  defect  problem.  Diaz-Bobillo’s  method  to  handle  this  issue  is  to 
replace  the  original  discrete  D  matrix  with  the  continuous  D  matrix.  However,  this  results 
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in  a  new  discrete  system  which  can  be  significantly  different  than  the  original  discrete 
system. 

This  problem  was  solved  by  simply  forcing  the  user  to  input  discrete  systems 
where  has  full  column  rank  and  has  full  row  rank.  Effectively,  these  requirements 
ensure  that  all  control  usage  is  penalized  and  that  no  measurements  are  perfect.  These 
conditions  are  also  standard  requirements  for  the  state- space  solution  of  H„  problems 

presented  in  [DGKF89].  With  this  change,  several  transformations  back  and  forth  from 
continuous  to  discrete  time  were  eliminated  in  Diaz-Bobillo’s  code.  This  proved  to  be  a 
major  overhaul  of  the  computer  code,  however,  since  many  parts  of  the  old  routine 
depended  on  the  continuous  time  version  of  the  original  discrete  system.  Among  other 
changes,  new  routines  had  to  be  developed  to  perform  Youla  parameterizations  and  stable 
projections  of  discrete  systems. 

Since  much  of  the  code  had  to  be  rewritten  to  incorporate  the  above  changes, 
other  simple  modifications  were  also  included  to  make  the  routine  more  user-friendly  and 
efficient.  As  a  result,  the  main  routine  is  now  a  function  file  instead  of  a  series  of  scripts, 
which  frees  up  additional  memory  in  the  MATLAB'^'^  workspace  for  the  user.  Data 
structures  used  in  the  |i-Analysis  and  Synthesis  Toolbox,  which  are  more  efficient  and 
easier  to  use,  were  also  incorporated.  All  the  modifications  to  the  original  software  were 
sent  back  to  researchers  at  the  Massachusetts  Institute  of  Technology  (MIT)  who  are 
preparing  the  optimization  software  for  commercial  release. 

3.3.  Understanding  Optimization 

In  order  to  answer  the  question  of  how  best  to  use  1^  optimization,  it  is  important  to 
first  understand  what  optimization  is  trying  to  accomplish.  The  formal  mathematical 

definition  given  in  Chapter  2  is  not  important  here;  rather,  a  simple  conceptual  idea  of  how 
the  method  works  is  sufficient.  By  definition,  optimization  attempts  to  minimize  the 

absolute  sum  of  a  system’s  sampled  pulse  response.  Conceptually,  this  optimization 
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method  works  by  pushing  down  on  the  pulse  response  from  all  sides.  In  other  words, 
both  peak-to-peak  gains  and  long  pulse  responses  are  penalized  since  both  tend  to  increase 
the  absolute  sum. 

Since  the  primary  interest  in  this  chapter  is  how  best  to  use  fj  optimization  for 
tracking  problems,  it  is  instructive  to  examine  the  unit  pulse  and  step  responses  of  a  simple 
discrete  system.  Consider  the  continuous  system. 


(3.1) 


(3.2) 


No.  of  Samples 

Figure  3.1.  Pulse  response  of  a  discrete  system 
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The  one-norm  of  this  system  can  be  calculated  by  inspection.  The  total  sum  of  samples  1- 
4  in  Figure  3.1  appears  to  be  approximately  1.  Indeed,  the  one-norm  for  this  system  is  1. 
The  unit  step  response  of  the  system  in  (3.2)  is  shown  in  Figure  3.2. 
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Figure  3.2.  Step  response  of  a  discrete  system 

Notice  the  distinct  relationship  between  the  unit  pulse  response  and  the  unit  step  response. 
If  r  is  the  sampled  unit  step  response  and  h  is  the  sampled  unit  pulse  response,  then 

r(kT)  =  ;£h(j)  for  k  =  1,2,...,  (3.3) 

j=i 

where  T  is  the  sample  period.  This  relationship  implies  that  the  faster  the  pulse  response 
decays  to  zero,  the  quicker  the  step  response  reaches  its  steady-state  value.  Since 
optimization  penalizes  long  pulse  responses,  it  logically  also  penalizes  slow  unit  step 
responses.  This  fact  and  the  general  relationship  between  the  unit  pulse  and  unit  step 
responses  are  particularly  important  in  using  fj  optimization  for  tracking  problems. 
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3.4.  Weighted  and  Unweighted  Sensitivity  Minimization 

The  goal  of  most  tracking  problems  is  to  minimize  the  error  between  the 
commanded  input  and  the  system  output.  This  type  of  problem  can  be  posed  as  a 
sensitivity  minimization  problem,  such  as  the  one  depicted  in  Figure  3.3. 

m 


Figure  3.3.  sensitivity  block  diagram 


For  the  AFTI  F-16  problem,  the  exogenous  input,  r,  is  an  unknown  commanded  normal 
acceleration  input  with  maximum  magnitude  less  than  or  equal  to  one,  and  the  controlled 
output,  m,  is  the  weighted  error  between  the  commanded  acceleration  and  the  actual 
aircraft  acceleration.  K  is  the  unknown  compensator,  and  a  state- space  description  of  the 
unweighted  plant,  G,  is  given  in  Appendix  A. 

For  now,  is  set  to  1  and  optimization  is  performed  on  the  above  system. 
Since  a  weighted  sensitivity  problem  places  no  penalty  on  control  usage,  a  small  penalty, 
p,  =  le  -  5 ,  is  added  to  to  ensure  the  left  inverse  of  U  exists.  The  optimal  closed-loop 
system  has  a  one-norm  of  2.07  and  the  controller  is  4th  order.  The  system  response  to  a 
Ig  step  in  normal  acceleration  is  shown  in  Figure  3.4.  The  “jags”  in  the  response  are  a 
product  of  the  sampled-data  simulation.  As  the  sample  rate  increases  these  “jags”  become 
less  apparent. 

Notice  that  the  step  response  is  extremely  fast.  This  is  mainly  due  to  the  fact  that 
there  was  only  a  small  penalty  placed  on  control  usage.  However,  as  discussed  in  the 
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Stabilator  deflection  (degrees)  Normal  acceleration  (g) 


previous  section,  unconstrained  l-^  optimization  tends  to  produce  very  quick  step 
responses.  Plots  of  control  usage  and  rate  of  control  usage  are  shown  in  Figures  3.5  and 
3.6. 


Time  (seconds) 

Figure  3.4.  Unweighted  sensitivity  step  response 


Figure  3.5.  Unweighted  sensitivity  control  usage 
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Figure  3.6.  Unweighted  sensitivity  control  rate 


The  control  usage  does  not  violate  the  maximum  deflection  limits,  but  it  is  quite  large  for 
only  a  Ig  change  in  normal  acceleration.  Since  the  system  is  linear,  it  is  easy  to  see  that 
the  maximum  deflection  limit  would  be  violated  for  a  commanded  2g  change  in  normal 
acceleration.  The  control  rate  violates  the  maximum  allowable  rate  limitation,  even  for  a 
very  small  command.  The  controller  found  above  would  be  undesirable  for  two  reasons: 
first,  the  system  tracks  with  a  steady-state  error;  second,  the  level  of  performance  shown 
in  Figure  3.4  is  unattainable  by  the  AFTI  F-16  due  to  limitations  in  the  stabilator  rate  of 
deflection. 

The  first  problem  can  be  handled  in  two  different  ways.  One  option  is  to  place  a 
gain  on  the  input  to  the  system  to  ensure  the  closed-loop  DC  gain  equals  one.  Another 
alternative  is  to  weight  the  sensitivity  more  heavily  at  low  frequency.  This  method  is 
preferred  since  it  can  eliminate  steady-state  error  problems  to  additional  low  frequency 
commands. 

To  demonstrate  the  second  method,  consider  the  following  weighting: 


W  = 


s-HO 
s-h  0.000 r 


(3.4) 
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This  particular  weighting  represents  the  inverse  of  the  desired  sensitivity  function  and  is 
very  similar  to  the  one  used  by  Luke  ([Luk93])  in  an  based  design  for  the  AFTI F-16. 

The  optimal  norm  for  the  system  using  this  weighting  is  2.25  and  the  compensator  is  5th 
order.  A  plot  of  the  system  response  to  a  Ig  step  in  normal  acceleration  is  shown  in 
Figure  3.7. 


Figure  3.7.  Weighted  sensitivity  step  response 


This  weighting  has  clearly  eliminated  the  problem  of  a  steady-state  error  to  a  step  input, 
but  the  system  response  is  still  extremely  fast.  Examination  of  the  control  usage  and  rate 
of  control  usage  for  this  weighting,  shown  in  Figures  3.8  and  3.9,  reveals  that  the  problem 
with  the  stabilator  deflection  rate  still  exists.  In  fact,  this  problem  is  worse  than  before. 

Before  discussing  how  to  handle  the  problem  of  excessive  control  deflection  and 
rates  of  deflection,  it  is  important  to  discuss  some  objectives  the  tracking  solution  should 
achieve.  The  following  list  represents  some  typical  factors  which  may  be  important: 

i)  minimum  error  to  low  frequency  commands 

ii)  no  violations  of  control  defection  and  rate  limitations 
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Stabilator  deflection  rate  (degrees/second)  Stabilator  deflection  (degrees) 


15 


Figure  3.8.  Weighted  sensitivity  control  usage 


Figure  3.9.  Weighted  sensitivity  control  rate 
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iii)  zero  steady-state  error  to  low  frequency  commands 

iv)  minimum  overshoot  and  undershoot 

v)  the  quickest  response  possible  given  the  above. 

Item  i)  indicates  that  sensitivity  minimization  is  the  proper  objective  function  for 
optimization,  but  items  iij-iv)  indicate  that  it  must  be  done  with  certain  constraints.  Item 
v)  is  inherently  built  into  optimization  for  most  problems.  To  demonstrate  how  items  i)- 
v)  can  be  achieved  without  frequency  weights  on  the  sensitivity  function,  is  set  equal  to 
one  in  the  following  sections. 

The  next  section  tackles  item  ii).  It  covers  three  different  approaches  for  adding 
control  deflection  and  rate  constraints  to  the  error  minimization  problem. 


3.5.  Control  Deflection  and  Rate  Limitations 

The  previous  section  was  concerned  with  solving  the  one-block  problem 


inf 

K  stabilizing 


(3.5) 


where  S  is  the  sensitivity  function.  This  section  is  concerned  with  the  general  two-block 
problem 


inf 

K  stabilizing 

where  W^.  is  a  weighting  on  the  control  usage.  The  added  block  in  (3.6)  can  be  used  to 
ensure  that  control  deflection  or  rate  limitations  are  not  violated. 

Since  the  control  rate  hmitations  were  violated  in  the  last  section,  only  rate 
constraints  are  added  in  this  section.  It  turns  out  that  once  the  control  rate  is  properly 
constrained  for  the  AFTI  F-16,  the  control  deflection  limitations  are  not  a  problem.  The 
ideas  presented  below  easily  extend  to  penalizing  control  deflections  alone  or  to  both 
control  rates  and  deflections. 


S 

W.KSI 


(3.6) 
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In  order  to  change  the  second  block  of  (3.6)  to  a  penalty  on  control  rate  instead  of 
control  usage,  an  appropriate  weight  must  be  chosen  for  W^.  Clearly  the  weighting  must 
be  chosen  so  that  it  effectively  takes  the  derivative  of  the  control  signal.  This  problem  is 
best  handled  directly  in  the  z-domain  with 

W,(2)  =  ^,  (3.7) 

Tz 

where  T  is  the  sample  period.  This  weighting  function,  known  as  the  backward  Euler 
transformation,  calculates  a  finite  different  gradient  between  discrete  pulses.  Since  the 
weighting  is  in  the  z  domain,  the  continuous  system  must  be  discretized  before  this  weight 
can  be  augmented  to  the  problem. 

The  first  approach  to  solving  the  rate-constrained  tracking  problem  is  to  multiply 
each  block  in  the  two  block  problem  by  a  desired  level  of  performance.  For  example,  if 
the  one-norm  of  the  first  block  is  desired  to  be  less  than  7  and  the  maximum  control 
deflection  rate  is  U^niax’  ^hen  the  problem  becomes 


inf 

K  stabilizing 


is 


u. 


-w  KS 


(3.8) 


If  the  resulting  one-norm  of  this  system  is  less  than  1,  then  both  levels  of  performance  are 
achieved.  This  follows  from  the  previous  assumption  that  the  maximum  magnitude  of  the 
exogenous  input  is  less  than  or  equal  to  one.  If  the  goal  is  to  find  a  solution  which  has  the 
minimum  achievable  7  without  violating  the  maximum  control  rate,  this  is  not  the  best 
approach  because  (3.8)  would  have  to  be  solved  iteratively  for  7  until  the  resulting  system 
one-norm  is  exactly  one. 

A  better  approach  is  to  solve  the  following  problem: 


inf  11S||^ 

K  stabilizing  " 

subject  to  ||W.KS||,<U„„. 
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Recall  for  multi-block  problems  that  the  one-norm  of  the  system  is  the  maximum  absolute 
row  sum.  In  (3.8),  the  maximum  absolute  row  sum  had  to  be  less  than  one  to  ensure  the 
one  norm  of  the  entire  system  was  also  less  than  one.  In  (3.9),  the  individual  row  sums 
are  separated,  with  one  being  minimized  while  the  other  is  constrained.  The  current 
version  of  the  software  developed  by  the  author  allows  the  user  to  specify  different 
constraint  levels  for  each  controlled  output  in  a  multi-block  column  problem. 

Problem  (3.9)  was  solved  for  the  AFTI  F-16  with  control  rate  limitations  of  200 
deg/sec,  100  deg/sec,  and  60  deg/sec.  A  plot  of  the  system  unit  step  response  for  all  three 
constraint  levels  is  shown  in  Figure  3.10. 


Figure  3.10.  Unweighted  sensitivity  step  responses  with  constraints  on  the  control  rate 

The  slowest  response  with  the  largest  steady-state  error  corresponds  to  the  actual 
stabilator  deflection  rate  limitation  of  60  deg/sec.  The  responses  to  the  other  two 
constraint  levels  are  shown  for  comparison.  This  type  of  plot  can  also  be  used  for  design 
purposes  since  it  is  easy  for  the  designer  to  see  how  much  performance  can  be  gained  if 
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faster  control  actuators  are  obtained.  Plots  of  the  control  deflections  and  rates  are  shown 


in  Figures  3.11  and  3.12,  respectively. 


Figure  3.11.  Unweighted  sensitivity  control  usage  with  constraints  on  the  control  rate 


Figure  3.12.  Unweighted  sensitivity  control  rates  with  fj  constraints  on  the  control  rate 
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Notice  from  Figure  3.12  that  the  control  rates  are  weU  within  their  respective  constraints 
for  a  unit  step  input.  The  constraints  will  actually  ensure  that  the  control  rate  limitations 

will  not  be  violated  for  any  input  with  a  maximum  magnitude  less  than  or  equal  to  one.  In 
this  particular  example,  the  constraints  are  very  conservative  for  a  unit  step  input. 

The  one-norm  of  the  objective  and  the  compensator  orders  are  shown  in  Table  3.1 
for  each  constraint  level. 


constraint 

one-norm 

controller 

order 

200  deg/sec 

2.63 

35 

100  deg/sec 

2.88 

38 

60  deg/sec 

3.17 

43 

Table  3.1.  Comparison  of  different  constraints  on  control  rate 

The  order  of  the  optimal  compensators  is  directly  related  to  the  support  length  of  the 
pulse  response;  i.e.,  the  number  of  time  steps  it  takes  the  pulse  response  to  decay  to  zero. 
In  fact,  if  the  optimization  software  did  not  perform  a  balanced  model  reduction  on  the 
compensator,  this  relationship  would  be  almost  one-to-one.  Recall  that  the  support  length 
of  the  pulse  response  is  related  to  the  time  it  takes  the  step  response  to  reach  steady-state. 
This  explains  why  the  controllers  found  using  the  above  approach  have  such  high  orders. 

The  previous  two  approaches  imposed  constraints  on  the  control  rate.  This 
means  that  the  constraint  limitation  imposed  will  not  be  exceeded  for  any  input  into  the 
system  bounded  in  magnitude  by  one.  Another  possibility  is  to  ensure  that  the  constraint 
is  not  exceeded  for  a  single  class  of  input  like  the  step  command.  Unlike  the  constraints, 
these  U  types  of  constraints  can  only  be  used  on  a  finite  horizon.  In  other  words,  they  can 
only  be  imposed  over  the  support  length  of  the  solution.  In  many  cases,  however, 
imposing  these  constraints  over  the  first  few  time  steps  is  sufficient. 
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To  understand  how  a  constraint  on  the  step  response  of  the  system  can  be  imposed 
in  optimization,  recall  the  relationship  in  (3.3)  between  the  unit  pulse  and  unit  step 

response.  The  step  response  at  any  particular  time  step  is  nothing  more  than  the  sum  of 
the  pulse  response  at  that  time  step  plus  all  previous  time  steps.  Therefore,  in  terms  of  the 
pulse  response  at  each  time  step,  these  constraints  can  be  imposed  with  very  simple 
Toeplitz  matrices,  with  ones  below  the  main  diagonal  and  zeros  above,  such  as 


- \ 

o 

o, 

'Mo)' 

1  0 

4>(i) 

< 

1 

:  ■•.  0 

• 

1 

1  1  1  i_ 

_6(N)_ 

i_ 

(3.10) 


Notice  that  these  constraints  can  easily  be  augmented  to  the  interpolation  conditions  in 
(2.41). 

The  new  problem  with  the  augmented  step  input  constraints  becomes 


inf  |lS|| 

K  Stabilizing  '  (3  11) 

subject  to  ||W„KSWf|L 

where  Wj  is  a  unit  step.  A  plot  of  the  AFTI  F-16  normal  acceleration  step  response  for 
control  rate  constraint  levels  of  200  deg/sec,  100  deg/sec  and  60  deg/sec  is  shown  in 
Figure  3.13. 

Again,  the  slowest  response  with  the  largest  steady-state  error  corresponds  to  a 
constraint  of  60  deg/sec.  Notice  that  the  step  responses  to  this  type  of  constraint  are 
much  quicker  and  have  less  steady-state  error  than  the  constraints.  Control  deflections 
are  shown  in  Figure  3.14  and  control  rates  are  shown  in  Figure  3.15. 

The  one-norm  of  the  objective  and  the  compensator  orders  are  shown  in  Table  3.2 
for  each  constraint  level.  With  this  approach,  quicker  settling  times  also  lead  to  lower 
order  controllers. 
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Figure  3.13.  Unweighted  sensitivity  step  responses  with  U  constraints  on  the  control  rate 


Figure  3.14.  Unweighted  sensitivity  control  usage  with  constraints  on  the  control  rate 
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Figure  3. 15.  Unweighted  sensitivity  control  rates  with  constraints  on  the  control  rate 


constraint 

one-norm 

controller 

order 

Mm 

2.16 

9 

100  deg/sec 

2.28 

13 

60  deg/sec 

2.43 

19 

Table  3.2.  Comparison  of  different  U  constraints  on  control  rate 

While  all  the  step  responses  in  this  section  meet  some  constraint  level  on  the 
control  rate,  none  of  them  has  zero  steady-state  error.  This  issue  is  addressed  in  the  next 
section. 

3.6.  Steady-State  Error  and  Time-Varying  Exponential  Weights 

Zero  steady-state  error  to  a  step  input  can  be  enforced  with  an  added  equality 
constraint  to  the  optimization  problem.  Recall  from  (3.3)  that  the  final  value  of  the  step 

response  is  simply  the  summation  of  the  unit  pulse  response  over  its  entire  support  length. 
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Therefore,  zero  steady-state  error  to  a  unit  step  input  is  guaranteed  if  the  sum  of  the 
sampled  unit  pulse  response  equals  zero.  This  is  not  an  absolute  summation  like  the  norm 
calculation;  it  is  simply  a  summation  of  the  pulse  response  at  each  time  step.  The  added 
equality  constraint  takes  the  form 


6(0) 

6(1) 


=  0. 


6(N) 


(3.12) 


This  constraint  was  added  to  the  problem  presented  in  (3. 11),  with  the  control  rate 
constraint  equal  to  60  deg/sec.  The  resulting  solution  has  an  objective  one-norm  of  3.00 
and  the  compensator  is  44th  order.  The  system  response  to  a  Ig  normal  acceleration  step 
input  is  shown  in  Figure  3. 16  and  the  control  rate  is  shown  in  Figure  3.17. 


Figure  3.16.  Unweighted  sensitivity  step  response  with  steady-state  error  and  control 
rate  constraints 
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Figure  3.17.  Unweighted  sensitivity  control  rate  with  steady-state  error  and  control  rate 
constraints 

Clearly,  zero  steady-state  error  is  achieved,  but  the  limit  on  control  rate  is  not  met  at  the 
nearly  discontinuous  jump  just  after  two  seconds.  Notice  that  the  response  in  Figure  3.16 
is  exactly  the  same  as  its  counterpart  in  Figure  3.13  up  until  this  jump.  This  happened 
because  the  steady-state  error  equality  constraint  was  not  imposed  until  the  very  last  time 
step  in  the  support  length.  In  this  system,  imposing  the  constraint  any  earlier  results  in  a 
larger  one-norm,  which  the  optimization  rejects. 

This  problem  can  be  overcome  with  time-varying  exponential  weights  on  the  norm 
calculation.  Consider  multiplying  each  sampled  pulse  response  by  a^”^,  where  k  is  the 
sample  index,  T  is  the  sample  period  and  a  >  1 .  Since  this  weight  gets  larger  as  k  gets 
larger,  it  effectively  penalizes  late  errors  over  early  ones. 

The  system  above  was  rerun  with  a  time-varying  exponential  weight  added  to  the 
norm  calculation.  With  a  =  1.1 ,  the  objective  one-norm  was  3.59  and  the  controller  was 
26th  order.  A  plot  of  the  system  response  to  a  step  input  is  shown  in  Figure  3.18  and  the 
control  rate  is  shown  in  Figure  3.19. 
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Figure  3.18.  Unweighted  sensitivity  step  response  with  steady-state  error/?^  control  rate 
constraints  and  time-varying  exponential  weights 
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Figure  3.19.  Unweighted  sensitivity  control  rate  with  steady-state  errorA^  control  rate 
constraints  and  time- varying  exponential  weights 
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Adding  the  exponential  weights  increased  the  one-norm  as  expected.  The  weighting  also 
decreased  the  settling  time  and  thus  the  controller  order. 

The  step  response  at  this  point  now  meets  all  the  criteria  established  in  Section  3.4, 
except  the  overshoot  issue.  This  problem  is  discussed  in  the  next  section. 

3.7.  Overshoot  and  Undershoot  Limitations 

Problems  with  excessive  overshoot  and  undershoot  in  the  step  response  can  be 
handled  in  exactly  the  same  manner  as  excessive  control  deflections  and  rate  violations. 
To  demonstrate  this  concept,  a  very  small  reduction  is  done  on  the  overshoot  for  the 
system  response  shown  in  Figure  3.18.  Theoretically,  the  overshoot  can  be  reduced  to  any 
desired  level  at  the  expense  of  a  slower  response.  However,  it  was  extremely  difficult  to 
find  a  solution  for  the  problem  presented  with  the  current  fj  optimization  software.  The 
multi-block  problem  contains  so  many  constraints  and  delays  that  calculation  attempts 
alone  are  extremely  expensive  in  terms  of  computer  time.  Further,  the  linear  programming 
routine  in  the  software  has  difficulty  solving  very  large  systems  of  equations,  possibly  due 
to  scaling  problems.  Both  of  these  issues  are  currently  being  addressed  by  researchers  at 
MIT. 

The  system  response  in  Figure  3.18  has  an  overshoot  of  about  80%.  An  type 
constraint  on  the  overshoot  was  added  to  the  problem  to  ensure  that  the  overshoot  would 
be  less  than  70%  to  a  step  input.  The  resulting  solution  had  an  objective  one-norm  of  3.40 
and  the  controller  was  28th  order.  A  plot  of  the  step  response  with  the  added  constraint  is 
shown  in  Figure  3.20. 

This  step  response  is  still  less  than  ideal;  however,  all  the  tools  to  shape  and 
constrain  the  response  are  now  available.  As  the  optimization  software  becomes  more 
reliable  and  efficient,  a  designer  should  be  able  to  use  all  the  techniques  presented  up  to 
this  point  to  find  a  compensator  which  meets  all  of  his  or  her  tracking  requirements. 
Unfortunately,  this  compensator  will  probably  be  extremely  high  order  and  therefore 
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impractical  to  use.  The  next  section  on  model  matching  offers  one  way  to  use 
optimization  and  still  produce  controllers  at  or  about  the  order  of  the  original  discrete 
system. 


Figure  3.20.  Unweighted  sensitivity  with  control  rate/steady- state  error/overshoot 
constraints  and  time-varying  exponential  weights 

3.8.  Model  Matching 

Using  optimization  to  solve  multi-block  systems  generally  produces  high  order 
controllers,  especially  when  additional  constraints  are  imposed  on  the  system.  One  way  to 
counter  this  problem  is  to  solve  a  one-block  system  instead.  Since  these  systems  can  be 
solved  exactly,  without  delay  augmentation,  the  resulting  controllers  tend  to  be  much 
smaller  (usually  about  the  order  of  the  unweighted  plant).  An  added  benefit  to  using  one- 
block  systems  is  that  they  can  be  solved  much  faster  and  more  reliably  than  multi-block 
systems  with  the  current  optimization  software. 

In  many  cases,  the  constraints  discussed  in  the  previous  sections  can  not  be 
imposed  with  one-block  systems.  Therefore,  a  one-block  system  must  be  chosen  that 
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incorporates  as  many  of  the  design  criteria  required  for  good  tracking  as  possible.  One 
way  to  accomplish  this  objective  is  to  model  match  the  design  problem  to  a  system  which 
has  the  desired  tracking  characteristics. 

Two  examples  of  model-matching  designs  are  presented  in  this  section.  The  first 
example  involves  the  SISO  AFTI  F-16  problem  that  has  been  discussed  throughout  this 
chapter.  The  second  example  involves  a  MIMO  missile  problem.  A  small  penalty  on 
control  usage,  similar  to  the  one  discussed  in  Section  3.4,  was  added  to  each  example  to 
make  the  resulting  problem  nonsingular. 

3.8.1.  SISO  Example 

Consider  the  closed-loop  model  matching  system  shown  in  Figure  3.21,  where  H  is 
the  ideal  closed-loop  model  given  in  continuous  time  by 

H(s)  =  -^  (3.13) 

s-(-4 


Figure  3.21.  Closed-loop  model  matching  diagram 


This  particular  closed-loop  model  was  chosen  because  its  step  response  is  relatively  quick, 
has  no  overshoot,  and  no  steady-state  error.  A  plot  of  the  step  response  of  this  model  is 
shown  in  Figure  3.22. 
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The  one-norm  of  the  solution  to  this  design  problem  is  0.37  and  the  compensator 
is  5th  order.  A  plot  of  the  AFTI  F-16  step  response  to  a  commanded  Ig  normal 
acceleration  change  is  shown  in  Figure  3.23. 


Figure  3.22.  Ideal  model  step  response 


Figure  3.23.  SISO  model  matching,  step  response 


3-25 


The  step  response  of  the  solution  is  approximately  0.37  g’s  larger  than  the  step  response 
of  the  ideal  model  at  steady-state.  To  make  the  system  response  match  the  ideal  model’s 
response,  the  commanded  normal  acceleration  has  to  be  multiplied  by  a  gain.  This  gain  is 
the  reciprocal  of  the  DC  gain  of  the  discrete  closed-loop  transfer  function, 


T^(z) 


K(z)G(z) 

H-K(z)G(z)‘ 


(3.14) 


For  this  problem,  the  gain  equals  0.73. 

Notice  that  with  this  particular  approach  there  is  no  direct  way  to  ensure  the  above 
solution  will  not  violate  control  deflection  and  rate  limitations.  In  this  problem,  the 
control  rate  limits  were  violated  in  the  first  few  time  steps.  However,  the  closed-loop 
system  still  performs  well  if  the  step  input  is  first  passed  through  a  continuous  prefilter 
equal  to 

F(s)  =  -^,  (3.15) 

s-i-10 

and  a  rate  limiter  is  added  to  the  control  signal.  With  the  added  prefilter,  the  system 
actually  sees  an  input  which  is  similar  to  the  response  in  Figure  3.22  (just  slightly  faster), 
instead  of  a  discontinuous  step  input.  The  new  input  is  actually  a  more  realistic 
representation  of  a  pilot  command. 

The  system  was  tested  with  the  prefilter,  gain  adjustments  on  the  input,  and  a 
control  rate  limiter  set  at  ±60  deg/sec.  The  response  is  shown  in  Figure  3.24.  This 
response  has  no  overshoot,  no  steady-state  error  and  was  achieved  with  a  5th  order 
controller  and  a  small  gain  on  the  input. 

Plots  of  the  sensitivity  and  complementary  sensitivity  for  the  resulting  closed-loop 
system  are  shown  in  Figures  3.25  and  3.26.  The  vector  gain  margins  (VGM)  and  phase 
margins  (VPM)  of  this  system  are 

-1 1.4  dB  <  VGM  <  10.4  dB  VPM  =  ±42.9° . 
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Figure  3.24.  SISO  model  matching,  step  response  with  scaling  factor,  prefilter  and 


control  rate  limiter 


Figure  3.25.  SISO  model  matching,  sensitivity 
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Figure  3.26.  SISO  model  matching,  complementary  sensitivity 

3.8.2.  MIMO  Example 

This  model  matching  example  is  done  on  the  3  input,  3  output  missile  system 
described  in  Appendix  B.  The  goal  of  this  design  problem  is  to  command  normal  or 
lateral  acceleration  while  keeping  the  roll-rate  of  the  missile  as  small  as  possible.  When 
commanding  normal  acceleration,  it  is  also  desired  to  minimize  the  lateral  acceleration  and 
vice  versa  when  commanding  lateral  acceleration. 

An  open-loop  model  matching  design  is  used  to  solve  this  problem,  as  shown  in 
Figure  3.27. 


Figure  3.27.  Open-loop  model  matching  diagram 
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In  this  problem,  the  ideal  open-loop  model  is  given  in  continuous  time  by 


was  chosen  to  determine  the  speed  of  the  responses.  Poles  were  placed  at  -0.001  in  each 
transfer  function  to  provide  near  integral  action,  and  ensure  very  little  steady-state  error. 

The  one-norm  of  the  solution  to  this  problem  is  1.009  and  the  controller  is  11th 
order.  To  avoid  control  rate  violations,  the  system  was  tested  with  rate  limiters  set  at  the 
maximum  allowable  rate,  750  degrees/second.  The  control  deflection  limits  (38  degrees) 


were  not  violated  in  any  case.  The  system  responses  to  a  lOg  normal  acceleration  step 
input  (applied  at  0.5  seconds)  are  shown  in  Figures  3.30  and  3.31.  The  system  responses 
to  a  5g  lateral  acceleration  step  input  are  shown  in  Figures  3.32  and  3.33.  Closed-loop 
sensitivity  and  complementary  sensitivity  plots  are  shown  in  Figures  3.28  and  3.29, 
respectively.  The  vector  gain  margins  (VGM)  and  phase  margins  (VPM)  of  this  system 


are 


-74.0  dB  <  VGM  <  1 1.2  dB  VPM  =  ±60.0° . 


Notice  that  the  open-loop  model  match  in  this  example  eliminates  the  need  for  the 
DC  gain  adjustment  used  in  the  SISO  example.  This  approach,  however,  lead  to  an 
unstable  system  in  the  AFTI F-16  problem  once  rate  limiters  were  added  to  the  simulation 
model.  While  neither  model  matching  method  is  guaranteed  to  produce  an  acceptable 
controller,  both  offer  easy  design  alternatives  to  the  methods  presented  in  Sections  3.4- 
3.7. 
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Figure  3.30.  MIMO  model  matching,  acceleration  responses  to  a  lOg  step  in  normal 
acceleration 


Figure  3.31.  MIMO  model  matching,  roll-rate  response  to  a  lOg  step  in  normal 
acceleration 
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Figure  3.32.  MIMO  model  matching,  acceleration  responses  to  a  5g  step  in  lateral 
acceleration 


Figure  3.33.  MIMO  model  matching,  roll-rate  response  to  5g  step  in  lateral  acceleration 
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Figure  3.34.  MIMO  model  matehing,  sensitivity 


Figure  3.35.  MIMO  model  matching,  complementary  sensitivity 
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3.9.  Chapter  Summary 

This  chapter  presented  methods  of  using  optimization  to  solve  tracking 
problems.  Several  constraints  that  could  be  added  to  the  standard  problem  to  handle 
control  deflection  and  rate  limitations,  zero  steady-state  error  requirements,  and  overshoot 
limitations  were  discussed.  Unfortunately,  the  constrained  optimization  problem  often 
produced  high  order  controllers.  To  counter  this  problem,  two  different  model  matching 
techniques  were  presented  which  can  produce  acceptable  tracking  results  with  lower  order 
controllers. 

The  next  design  issue  to  be  considered  is  system  noise  performance,  which  is  best 
handled  with  H2  optimization.  To  achieve  both  tracking  and  noise  performance,  one  might 
suggest  using  a  mix  of  H2  and  optimization  techniques.  Indeed,  this  is  a  viable  approach. 
Mixed  H2/LJ  optimization  for  continuous  time  systems  is  discuss  in  Chapter  4  and  mixed 
HjAi  optimization  for  discrete  time  systems  is  discussed  in  Chapter  5.  In  both  cases,  the 
order  of  the  resulting  compensator  is  chosen  by  the  designer  which  eliminates  the  problem 
of  optimal  solutions  with  excessively  large  compensator  orders. 
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N.  Mixed  HJL^  Optimization 


4.1.  Chapter  Overview 

This  chapter  develops  a  fixed-order,  mixed  optimization  method  to  design 
controllers  for  continuous  systems.  Section  4.2  discusses  the  general  mixed  H2/L1 
optimization  problem.  In  Section  4.3,  analytical  derivatives  for  the  portion  of  the 
problem  are  derived.  Two  methods  of  numerically  approximating  the  Lj  norm  of  a 
transfer  function  are  presented  in  Section  4.4.  Analytical  gradients  for  the  Lj  portion  of 
the  problem  and  a  stability  constraint  are  discussed  in  Section  4.5.  Section  4.6  covers  the 
computer  implementation  of  the  Hj/L,  problem.  The  last  section  in  the  chapter 
demonstrates  H^/Lj  optimization  using  a  continuous  version  of  the  AFTI F-16  longitudinal 
model. 


4.2.  The  Mixed  Problem 

The  mixed  H^/Lj  design  problem  is  depicted  in  Figure  4.1. 


W 

r 

u 


Figure  4. 1 .  Mixed  TIJL^  problem 


Z 

m 

y 


In  general,  it  is  assumed  that  there  is  no  relationship  between  w  and  r,  or  z  and  m. 

The  goal  of  mixed  optimization  is  to  find  a  stabilizing  controller  which 
satisfies 
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where 


(4.1) 


T,..,  = 


T  = 

mr 


Aa 

1 - 

C. 

DzJ’ 

A, 

Bri 

D^rj 

(4.2) 


(4.3) 


and  V  is  a  user  specified  constraint  level  on  the  Lj  norm  of  T^. 

The  state- space  of  P  is  found  by  augmenting  the  stable  weights  of  the  H2  problem 
and  the  Lj  problem  to  the  original  system.  Typically,  the  orders  of  the  individual  H2  and 
Lj  problems  are  less  than  the  order  of  P.  The  state-space  equations  of  the  Hj  and  Lj 
problems  can  be  written  as 


X2  =  A2X2  +B^W  +  B„2U 

z  =  C,X2+D,,w  +  D,„u 

y  =  Cy2X2+Dy„w-f-Dy„u 

Xi  =  AjXi  -HB,r  +  B„,u 
m=C„,Xi+D„„r-(-D^„u 

y  =  CyiXi-(-Dy,r-fD^u, 


(4.4) 


(4.5) 


where  X2  is  the  state  vector  for  the  underlying  H2  problem,  and  Xj  is  the  state  vector  of  the 
underlying  problem.  The  following  assumptions  are  now  made  on  the  state-space 
elements  in  (4.4)  and  (4.5): 

i)  Dzw=0 

ii)  D^=0 


iii)  (A2,  B^^)  is  stabilizable  and  (C^^,  is  detectable 


iv)  and  have  full  rank 


V) 


A2-j0)I  B„2 

C,  D,, 


has  full  column  rank  for  all  0) 
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vi) 


has  full  row  rank  for  all  co 


A2-jcoI  B„' 

C^2  Dy., 

vii)  (Ap  B„j)  is  stabilizable  and  (Cyp  Aj)  is  detectable 
Reasons  for  assumptions  i)-vi)  were  given  in  Section  2.8.  Assumption  vii)  is  required  to 
ensure  the  Lj  problem  has  a  solution. 

The  controller  state- space  equations  are 

u=Qx^-l-D^y.  (4.6) 

Using  (4.4)  and  (4.6),  the  closed-loop  state-space  equations  for  can  be  written  as 

Xj  =  (A2  +  B^2D^Cy2)x2  +  By2CjjX^  +  (B^  -f-  B^2Dj.Dy^)w 
Xk  =  BkCy2X2  +  A.Xj^  B.Dy^w 

z  =  (C,  -fD^DkCy2)x2  +  D^Qx,  -hD^Dj^Dy^w.  (4.7) 

Notice  that  D^D^Dy^  in  (4.7)  must  be  zero  to  ensure  the  resulting  two-norm  of  is 
finite.  This  fact  and  assumption  iv)  imply  that  Dj^  must  be  zero.  Therefore,  a  strictly 
proper  controller,  K,  can  be  assumed  without  any  loss  of  generality.  With  this  additional 

assumption,  the  state-space  equations  for  the  mixed  H^/L^  problem  can  be  written  as 

X2  =  A2X2-»-B„W 

Z —  Cj,X2 
Xj  =  AiXj-l-B.r 

m=C„Xj-J-D„,r,  (4.8) 

where 


(4.10) 


A2  ®u2^k 

®k^y2  Ajj 


(4.11) 
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^1  ^ul^k 

®k^yl  ^k 


(4.12) 


B 


w 


BkDyw 


(4.13) 


B,  = 


B, 

BkD, 


C.=[C,  D^CJ 
C^=K  D.„CJ 

D_  =D„.. 


(4.14) 

(4.15) 

(4.16) 

(4.17) 


The  following  definitions  are  used  to  discuss  the  solution  to  the  mixed  problem: 


V  =  inf  ||T„|I 

K  admissible*’  '  ^ 


K2  opt  =  the  unique  K  that  makes  IjT^^  =  a 
Ki  op  =  a  K  that  makes  ||T^[  =  v 
®  =  ItlL  when  K=K,, 


[  Opt 

V  =  ||T_JL  whenK=K, 

ll  mr||i  z  opt 


Kmp  =  the  global  minimum  solution  to  the  H2  /  Lj  problem  for  some  v 


«  ^  F^wi,  whenK  =  K^p 

v‘=  \Kl  whenK^K^P 


(4.18) 


Admissible  controllers  must  be  stabilizing  and  have  a  fixed,  user- specified  order. 

A  solution  to  the  YIJL^  problem  must  satisfy  the  Kuhn-Tucker  necessary 


conditions: 


i)  must  be  feasible,  i.e.  stabilize  and 
“)^tT„IL+)^.V||T4=0,  k,>0 
m)k,(||T4-v)=0,k,>0 
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where  Xj  is  a  Lagrange  multiplier  associated  with  the  one-norm  constraint.  Condition  i)  is 
simply  a  feasibility  condition.  Condition  ii)  states  that  the  gradient  of  the  objective  must 
be  balanced  by  the  scaled  gradient  of  the  constraint.  The  last  condition  states  that  if  the 
constraint  is  not  satisfied  exactly,  then  A,,  must  be  zero.  Conditions  ii)  and  iii)  together 

imply  that  an  optimal  solution  (if  it  exists)  must  lie  on  the  constraint  boundary  when 
v<v<v.  Ifv>v,  then  the  unique  solution  is  K2opt .  By  definition,  v  can  not  be  chosen 

less  than  v . 

When  V  <  V  <  V ,  a*  is  a  monotonically  decreasing  function  of  v.  An  illustration  of 
a  possible  Hj/Lj  solution  curve  is  shown  in  Figure  4.2. 


V 

Figure  4.2.  Mixed  solution  boundary 

A  sequential  quadratic  programming  (SQP)  algorithm  is  used  to  solve  the 
problem  numerically.  The  objective  (f)  and  the  constraints  (g)  are  given  by 

«c)  =  S.lT„|^ 

g.w=e.(|rr„||,-v) 

g.(c)  =  S.{m|K(Re(>-i(A2)))|,  (4.19) 
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where  c  is  a  vectorized  compensator,  q’s  are  scaling  parameters,  and  ?i;(  )  is  the  i* 
eigenvalue  of  (•) .  The  stability  constraint,  g^,  is  added  to  the  problem  to  keep  the  SQP 

algorithm  from  getting  lost  in  the  unstable  region.  While  this  constraint  should  be  posed 
as  a  strict  inequality  constraint,  equality  must  also  be  allowed  to  incorporate  it  into  the 
SQP  algorithm.  Modal  form  is  assumed  for  the  controller,  K,  to  minimize  the  number  of 
design  variables  in  c.  While  this  approach  disallows  repeated  eigenvalues  in  the  controller, 
it  has  been  shown  to  be  sufficient  in  practice. 

The  SQP  algorithm  requires  gradients  for  the  objective  function  and  each  of  the 
constraints.  Analytical  gradients  for  the  objective  function  are  derived  in  the  next  section. 

4.3.  Gradients  of  the  T wo  -Norm 

Recall  from  Section  2.8  that  the  two-norm  squared  of  T^^  is  given  by 

||T..||’  =  tr[QClC.],  (4.20) 

where  Q  is  the  positive  semidefinite  solution  to  the  Lyapunov  equation 

A2Q-f-QA^-KBX=0.  (4.21) 

If  the  constraint  on  Q  in  (4.21)  is  written  in  an  equivalent  form, 

tr[(A2Q4-  QA^  +  B„B^)X]  =  0  for  all  X,  (4.22) 

then 

||T^  1^  =  tr[QC^C  J tr[(A,Q+ QA^  +B„B:)X]  for  all  X .  (4.23) 

To  simplify  the  notation,  let  J  equal  the  two-norm  squared  of  T^,^.  Notice  from  (4.23)  that 
J  is  an  explicit  function  of  the  design  variables  and  the  matrices,  X  and  Q.  However,  the 
partial  derivative  of  J  with  respect  to  X  is  simply  the  left  side  of  (4.21)  and,  therefore,  it  is 
always  equal  to  zero.  Further,  it  can  be  shown  that 

^=A;x+XA,  +  C;C.,  (4.24) 
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which  is  a  Lyapunov  equation.  Since  X  is  arbitrary,  let  X  be  the  positive  semidefinite 
solution  to  the  Lyapunov  equation 

A^X  +  XA2+CJC,=0.  (4.25) 

This  choice  of  X  is  guaranteed  to  exist  since  Aj  must  be  stable.  With  this  choice,  the 
partial  of  J  with  respect  to  Q  is  also  zero,  and  the  only  remaining  derivatives  to  be 
calculated  are  the  partials  of  J  with  respect  to  the  design  variables  (i.e.  A,^,  and  CJ. 

If  Q  is  given  by 


and  X  is  given  by 


then  these  derivatives  are 


X  = 


X„  X 


12 


lXu  x,2, 


(4.26) 


(4.27) 


aj 

aj 

aBk 

aj 

ac. 


—  2[^Xi2Qi2  +  X22Q22] 

=  2[x^2QnCj2  +X22Q?2C^2  +X[2B,D;„  +X2B,D^DJ,] 
“  ^[®u2XiiQi2  +  By2Xi2Q22  +  DzuCzQl2  "*■  ^zu^zuCkQ22]- 


(4.28) 


The  complete  method  for  finding  the  partial  derivatives  of  J  with  respect  to  the 
vectorized  compensator,  c,  is  as  follows: 

i)  Solve  the  Lyapunov  equations  in  (4.21)  and  (4.25)  for  Q  and  X  respectively 

ii)  Compute  the  partial  derivatives  in  (4.28) 

iii)  Rearrange  (4.28)  to  express  the  gradient  with  respect  to  the  vector  of  design 
variables,  c,  as  a  vector. 
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£ 

9c 

where  the  individual  vectors  in  parentheses  are  the  columns  of  the  partial  derivative 
matrices,  %  is  the  order  of  the  compensator,  and  n^  is  the  number  of  measurements.  The 
next  section  discusses  two  numerical  methods  of  approximating  the  Lj  norm  of  T^. 

4.4.  Calculating  the  Lj  norm 

As  mentioned  in  Section  2.7,  there  is  no  known  way  of  analytically  computing  the 
Lj  norm  of  a  continuous  system.  Many  of  the  current  approaches  to  the  problem  involve 
discretizing  the  continuous  system  with  the  forward  Euler  rule,  and  computing  the  fj  norm 
of  the  resulting  discrete  system.  These  methods,  however,  prove  to  be  impractical  in 
many  cases.  This  section  presents  two  numerical  approaches  which  attack  the  problem 
directly,  and  avoid  the  complications  associated  with  forward  Euler  transformation. 

The  Lj  norm  of  a  continuous,  SISO  system  is  given  by 

oo 

l|Tj.=J|C„e‘-B,|dt+|D^|,  (4.30) 

0 

If  Aj  is  stable,  then  e'^'*  approaches  zero  as  t  approaches  infinity,  and  (4.30)  can  be 
approximated  by  tmncating  the  integral  at  some  point  in  time,  tj,,,  which  can  be  computed 
from  the  eigenvalues  of  A^.  The  remaining  question  is  how  best  to  evaluate  the  truncated 
integral. 

The  first  approach  to  this  problem  is  to  eliminate  the  absolute  value  sign  inside  the 
integral  by  determining  where  the  impulse  response  is  positive  and  negative.  This  can  be 
done  by  discretizing  the  continuous  system  with  a  ZOH,  and  finding  the  pulse  response  of 
the  resulting  discrete  system.  Since  the  ZOH  transformation  preserves  the  values  of  the 
continuous  impulse  response  at  the  sample  points,  the  discrete  pulse  response  can  be  used 


JL] 

Ji 


aj 


"k 


aj  ] 

v^bJj 


^  aj 


v^^k  y„ 


^aj 


ni 


(4.29) 


4-8 


to  find  approximate  locations  where  the  continuous  impulse  response  is  zero.  These 
approximate  locations  are  then  refined  to  any  degree  of  accuracy  by  using  a  nonlinear 
root-solver  on  the  continuous  function,  ^ .  Once  the  zero  locations  of  the  impulse 

response  are  known,  the  absolute  value  sign  can  be  removed  by  breaking  the  integral  in 
(4.30)  into  a  series  of  integrations.  Further,  if  Aj  is  invertible  (i.e.  has  no  zero 
eigenvalues),  then 

^2 

jC„e''‘^B,dt  =  C„A“'[e''*  -2I]B, .  (4.31) 

Thus,  each  integral  in  the  series  can  be  determined  exactly.  When  A^  does  have 
eigenvalues  at  the  origin,  the  system  is  unstable  and  the  one-norm  is  infinite. 

The  key  to  this  method  rests  on  determining  the  zero  locations  of  the  continuous 
impulse  response.  If  zero  locations  are  missed  in  the  discretization  step  due  to  high 
frequency  dynamics,  then  the  one-norm  wUl  be  inaccurate.  In  addition,  most  nonlinear 
root-solvers  are  only  capable  of  finding  the  closest  root  to  a  given  initial  guess.  This 
implies  that  the  approximate  root  locations  must  be  fairly  precise,  i.e.  the  discretization 
sample  period  must  be  fairly  small.  Using  a  small  sample  period  to  estimate  the  impulse 
response  of  systems  with  fast  and  slow  dynamics  can  be  expensive  in  terms  of  computer 
time.  These  issues  make  the  above  method  impractical  to  use  in  a  computer  algorithm 
which  must  handle  a  wide  range  of  problems.  However,  this  method  does  offer  an 
alternative  way  to  calculate  the  one-norm  for  specific  examples,  and  can  be  used  to  check 
another  method. 

The  second  approach  to  approximating  the  truncated  integral  is  more  robust  but 
less  accurate.  In  this  method,  the  continuous  system  is  discretized  with  a  ZOH 
transformation  using  a  small  sample  period.  The  truncated  integral  is  then  approximated 
with  a  trapezoidal  integration  of  the  discrete  pulse  response.  As  the  sample  period 
decreases,  the  approximation  improves.  If  H/L^  optimization  is  performed  near  Lj 
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optimal,  where  is  small,  this  method  works  well  at  approximating  the  one-norm  of  a 
system  without  requiring  an  unreasonable  number  of  samples.  However,  finding  a  starting 
point  near  Lj  optimal  can  be  difficult.  Currently,  the  most  reliable  method  of  finding  an 
adequate  starting  point  is  to  use  a  solution  from  fixed-order  H^/H^  optimization  ([Wal94]) 
near  optimal  (see  [Rid92]  for  an  excellent  discussion  of  H„  optimization). 

The  Lj  norm  for  MIMO  systems  is  calculated  from  the  maximum  row  sum  of  SISO 
transfer  function  norms.  However,  discontinuous  gradients  can  occur  when  the  maximum 
row  sum  occurs  over  more  than  one  row.  To  counter  this  problem,  each  row  sum  is 
constrained  as  a  separate  Multiple  Input  Single  Output  (MISO)  transfer  function.  This 
effectively  adds  more  constraints  to  the  optimization  problem,  but  most  of  the 
additional  constraints  are  not  active  at  any  specific  design  point.  If  an  optimization 
algorithm  is  used  which  only  calls  for  gradients  of  the  active  constraints,  then  these 
additional  constraints  have  little  impact  on  the  overall  performance  of  the  Hj/Lj  algorithm. 
The  next  section  discusses  how  to  calculate  the  gradients  associated  with  the  Lj  norm. 

4.5.  Gradients  of  the  One-Norm 

This  section  derives  the  gradients  of  the  one-norm  with  respect  to  the  design 
variables  in  the  matrices  Aj^,  B^,  and  Q.  Since  all  of  these  matrices  appear  in  Aj  [see 
(4.12)],  the  first  problem  considered  is  how  to  compute  the  partial  of  e^*‘  with  respect  to 
any  element  of  A^.  This  is  the  most  difficult  issue  in  calculating  the  gradients  of  the  one- 
norm  with  respect  to  the  design  variables.  Contributions  to  the  one-norm  gradients  from 
the  other  closed-loop  state-space  matrices,  which  are  considerably  simpler  to  find,  are 
discussed  thereafter.  In  solving  the  first  problem,  it  wiU  also  become  clear  how  to 
calculate  the  gradient  of  the  added  stability  constraint  in  (4.19). 

If  Aj  is  a  non-defective  matrix  (Aj  can  be  diagonalized),  then  the  partial  of  e'^^' 
with  respect  to  any  element  of  Aj  is  given  by 
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(4.32) 


i^  =  AfRe-R-]=®e"R-+R^R-‘  +  Re":^, 
3a-  3a  -  ^  ^  3a  3a  3a- 

ij  **ij  ^‘*ij 


da,.  da,. 


where  R  is  the  right  eigenvector  matrix  of  Aj  and  A  is  a  diagonal  matrix  consisting  of  the 


eigenvalues  of  A^.  The  partial  of  e^‘  with  respect  to  a^.  can  be  computed  easily  from 


ae^‘_Aae^‘  dX, 

da, 


(4.33) 


provided  that  the  partial  derivative  of  each  eigenvalue  is  known  with  respect  to  a,. 
Likewise,  the  partial  derivative  of  R  *  with  respect  to  a,  can  be  easily  computed  from 


aR-‘_  „-.3R„-, 


(4.34) 


provided  that  the  partial  derivative  of  the  right  eigenvector  matrix  with  respect  to  a,  is 
known. 

The  partials  of  an  eigenvalue  and  eigenvector  with  respect  to  ajj  can  be  computed 
from  the  standard  eigenvalue  problem  if  the  Uj  eigenvalues  of  A^  are  distinct.  The 
derivations  of  these  two  derivatives  involve  both  the  left  and  right  eigenvectors  of  A^,  and 
are  described  in  detail  in  [Nel76].  Letting 


the  solutions  are  given  by 


da, 


X'=l[a;r, 


(4.35) 


R'  -  +CiRi  -  X  H-CjRi, 


(4.36) 


where 


Ll[R,k'-A;Ri] 


,  i  ^  k 


(4.37) 
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(4.38) 


Cj  =  -Re(Rf  MVj-^Rf  MX 
and 

RfMR.  =  l.  (4.39) 

The  symbol,  (0^,  denotes  the  complex  conjugate  transpose  of  (•),  and  L  is  the  left 
eigenvector  matrix.  R^  and  L-  refer  to  the  i*  eigenvector  of  R  and  L,  respectively.  Note 
that  A[  is  simply  a  Uj  x  nj  matrix  of  zeros  with  a  one  in  the  (i,j)  element.  The  real  part  of 
(4.35)  provides  a  method  of  computing  the  gradients  of  the  stability  constraint  in  (4.19)  if 
is  the  maximum  eigenvalue  of  A^. 

The  partial  of  e'^*^  with  respect  to  3i-  can  be  completely  determined  from  (4.33)- 
(4.39).  Using  this  information,  the  partial  derivative  of  the  one-norm  with  respect  to  Aj 
can  be  found  element-by-element  from 


(4.40) 


where  sgn(0  is  1,  -1,  or  0  depending  on  the  sign  of  (•).  The  partial  derivatives  with 
respect  to  the  other  closed-loop  matrices  are  given  by 

3||T. 


mrlli 


3B. 


=  |sgn(C„e*-B,)(e*'‘)  C 


dt 


(4.41) 


Jlgn(C„e^^'B,)[B^(e^‘f 


ac„ 


dt 


(4.42) 


30. 


=  sgn(D„,). 


(4.43) 


From  these  expressions,  the  partials  of  the  one-norm  with  respect  to  Aj^,  Bj^,  and  Q  can  be 
expressed  as 
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fa_ 

3A 


4U 

0Aj 


(4.44) 


Jiij  +i,ni+j 


3B  y^M 


p=l  q=l 


au 

3Aj 


+  ££^ppyr. 


nj+i,n|+j 


p=l  q=l 


3||T.J, 


36. 


(4.45) 


Jrlj+p.q 


3||T. 


3C. 


^  =  SX5,,B,, 


P=1  q=l 


'3|iT. 


mrlli 


3Ai 


Jp,ni+q 


+  1X5, „D„, 

p=l  q=I 


5|tlL 


3C„ 


(4.46) 


Jp,n,+q 


The  partial  of  the  one-norm  of  with  respect  to  c  can  be  formed  from  the  columns  of 
matrices  in  (4.44)-(4.46)  as  described  in  Section  4.3. 

The  integrals  in  (4.40)-(4.43)  are  approximated  in  the  same  manner  as  the  integral 
in  the  one-norm  calculation.  The  upper  limit  of  integration  is  changed  to  tjq,  and 
trapezoidal  integration  is  performed  over  discrete  points  from  the  continuous  function. 
The  same  discretization  period  is  used  for  the  norm  calculation  and  the  gradients,  which 
allows  the  sign  factor  required  m  (4.40)-(4.43)  to  be  computed  from  the  data  gathered  in 
evaluating  the  one-norm. 

As  mentioned  in  the  previous  section,  separate  MISO  gradients  are  calculated  for 
each  row  of  a  MIMO  system.  This  ensures  that  continuous  gradient  information  is 
available  regardless  of  where  the  maximum  row  sum  occurs. 

Recall  that  development  of  R'  and  X'  require  that  the  eigenvalues  of  A^  be 
distinct.  Clearly  this  condition  can  be  violated  while  performing  optimization.  The 
current  algorithm  switches  to  finite  difference  calculations  for  gradient  information  if  this 
occurs.  Research  is  still  being  conducted  to  develop  ways  to  handle  an  A^  matrix  with 
repeated  eigenvalues. 


4-13 


4.6.  Computer  Implementation 

The  H2/L1  optimization  problem  is  currently  implemented  using  the  MATLAB'^'^ 
SQP  routine,  constr.m.  Separate  subroutines  for  each  norm,  constraint  and  gradient 
calculation  are  called  from  this  routine.  Unfortunately,  this  routine  currently  requires  all 
gradient  information  regardless  of  whether  or  not  a  constraint  is  active.  FORTRAhT” 
shells  can  be  written  to  allow  the  MATLAB'^“  subroutines  to  interface  with  IMSL^'^ 
optimization  routines,  which  eliminate  this  problem.  Section  4.7  demonstrates  the  H2/LJ 
optimization  algorithm  on  a  design  problem  for  the  AFTI F-16. 

4.7.  H^IL^  Design  Example 

A  longitudinal  controller  design  for  the  AFTI  F-16,  shown  in  detail  in  Appendix  A, 
is  used  to  demonstrate  the  H^/Lj  optimization  method.  The  primary  objective  of  the  Hj 
portion  of  the  problem  is  to  minimize  the  effect  of  wind  disturbances  and  measurement 
noise.  The  Lj  portion  of  the  problem  is  designed  to  improve  tracking  performance. 

4.7.1.  The  H2  Problem 

The  H2  problem  for  the  AFTI  F-16,  taken  from  [Luk93],  is  to  find  an  internally 
stabilizing  controller  which  minimizes  the  response  of  the  normal  acceleration  and 
weighted  control  to  the  wind  disturbances  and  measurement  noise.  A  block  diagram  of 
the  problem  is  shown  in  Figure  4.3. 
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Figure  4.3.  F- 1 6  H2  diagram 


The  control  weight,  p,  equals  10  and  the  state  weighting  matrix,  H,  equals  the 
system  C  matrix.  The  wind  disturbance,  Wj,  is  modeled  as  a  white  Gaussian  noise  (WGN) 
with  5.0  X 10^  rad^-sec  intensity,  and  F  is  the  column  of  A  corresponding  to  the  angle-of- 
attack  state,  a.  The  measurement  noise,  W2,  is  modeled  as  a  WGN  with  1.6  x  10'^  g^-see 
intensity  and  w^  equals  1.  Weighted  control  power,  Zj,  and  normal  acceleration,  z^,  are 
the  controlled  outputs.  The  plant  state-space  matrices  are  given  in  Appendix  A  and  the 
resulting  matrices  are  given  in  Appendix  C. 


4.7.2.  The  Lj  Problem 

The  L;  problem,  depicted  in  Figure  4.4,  is  a  weighted  sensitivity  minimization 
design.  Wj  is  the  inverse  of  the  desired  sensitivity,  and  is  given  by 

s+10 


W  =■ 


s+0.0001 


(4.47) 


A  plot  of  the  desired  sensitivity  is  shown  m  Figure  4.5.  The  matrices  for  the  Lj  problem 
are  also  given  in  Appendix  C. 
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m 


Figure  4.4.  F-16  diagram 


Figure  4.5.  F-16  desired  sensitivity 
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4.7.3.  Results 

A  plot  of  the  Hj/Lj  solution  curve  for  a  6*  order  compensator  is  shown  in  Figure 
4.6. 


Figure  4.6.  F-16  H2/L1  solution  curve 


Notice  that  the  plot  only  depicts  mixed  solutions  near  Lj  optimal.  A  complete 
solution  curve  would  extend  much  further  along  the  horizontal  axis,  and  asymptotically 
approach  the  dashed  line,  labeled  pc .  Values  for  several  of  the  solution  points,  depicted  as 

circles  in  Figure  4.6,  are  given  in  Table  4.1. 

The  points  are  numbered  from  left  to  right  in  Figure  4.6.  Point  number  24 
represents  the  H2  optimal  solution,  which  is  not  depicted  in  Figure  4.6.  The  points  in 
Table  4. 1  are  referred  to  as  the  design  points  for  the  remainder  of  this  section. 

Plots  of  the  sensitivity  and  complementary  sensitivity  for  the  design  points  are 
shown  in  Figures  4.7  and  4.8. 
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point  # 

a 

V 

1 

3.23 

3.40 

5 

2.40 

3.96 

9 

2.07 

4.36 

15 

1.67 

5.02 

21 

1.50 

5.66 

24 

0.86 

97,190 

Table  4.1.  F-16  solution  points 


Figure  4.7.  F-16  sensitivity 
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Frequency  (radians/second) 

Figure  4.8.  F-16  complementary  sensitivity 

In  each  plot,  one  curve  is  dramatically  different  than  the  others.  This  curve  corresponds  to 
design  point  24,  the  Flj  solution. 

The  vector  gain  and  phase  margins  for  the  design  points  are  shown  in  Table  4.2. 


point  # 

lower  VGM  (dB) 

upper  VGM  (dB) 

VPM  (deg) 

1 

-7.77 

5.42 

34.4 

5 

-10.3 

5.42 

40.7 

9 

-9.12 

6.10 

37.9 

15 

-7.88 

8.02 

35.1 

21 

-7.50 

8.11 

35.3 

24 

-7.82 

8.80 

37.1 

Table  4.2.  F-16  vector  gain  and  phase  margins 


Notice  that  the  stability  margins  do  not  consistently  improve  as  v  is  decreased.  However, 
the  margins  are  acceptable  at  all  of  the  design  points. 
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The  open-loop  GK  plot  is  shown  in  Figure  4.9. 


The  shaded  area  on  the  left  side  of  Figure  4.9  represents  a  performance  and 
disturbance  rejection  barrier.  The  shaded  area  on  the  right  side  of  Figure  4.9  represents  a 
sensor  noise  and  unmodeled  dynamics  barrier.  Descriptions  of  both  barriers  were  taken 
from  [RB86],  which  also  contains  an  excellent  discussion  of  desired  GK  shapes. 

Notice  that  the  mixed  design  points,  all  relatively  near  Lj  optimal,  miss  the  barrier 
on  the  left,  but  pass  within  the  barrier  on  the  right.  A  design  which  meets  both  barriers 
could  be  found  by  a  mixed  solution  closer  to  H2  optimal.  Alternatively,  the  mixed 
optimization  design  could  be  redone  with  the  weighting 

W3  = 

’  s-t-aoooi 

on  the  Lj  problem,  instead  of  the  weighting  in  (4.47).  This  new  weighting  should  provide 
acceptable  GK  loop  shapes  for  mixed  HL^/Lj  solutions  near  Lj  optimal,  without  decreasing 
the  bandwidth  of  the  system  below  the  frequency  of  average  pilot  commands. 
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The  loop  shapes  from  the  design  points  imply  that  high  frequency  noise  will  be 
more  prevalent  than  low  frequency  noise  in  system  responses.  Systems  obtained  from 
these  design  points  should  also  track  low  frequency  commands  well.  Both  of  these  facts 
can  be  seen  in  the  step  responses  presented  below. 

Figure  4.10  shows  the  F-16  responses,  without  noise,  to  a  commanded  Ig  step 
input  for  the  different  designs.  Notice  that  the  H2  solution  tracks  with  a  steady-state  error 
while  the  mixed  designs  do  not.  Step  responses  with  noise  for  design  points  1,  9,  15  and 
24  are  shown  in  Figures  4.11-4.14,  respectively.  As  expected,  the  systems  with  lower  v 
values  contain  more  high  frequency  noise  than  those  closer  to  the  H2  optimal  design.  The 
most  important  item  to  note  from  Figures  4.11-4.14  is  that  the  tracking  performance  of 
the  AFTI  F-16  can  be  greatly  improved  using  mixed  designs,  with  very  little  increase 
on  the  amount  of  noise  in  the  system  response. 


Figure  4.10.  F-16  step  responses  without  noise 
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Figure  4.1 1.  F-16  step  response  with  noise,  design  point  1 


Figure  4.12.  F-16  step  response  with  noise,  design  point  9 
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Figure  4.13.  F-16  step  response  with  noise,  design  point  15 


Figure  4. 14.  F-16  step  response  with  noise,  design  point  24 


Design  point  15  was  chosen  as  the  best  mixed  design  based  on  stability 
margins,  loop  shape  and  step  responses.  This  particular  design  will  be  compared  to  the 
best  discrete  design  in  the  next  chapter.  The  control  usage  for  this  particular  design 

point  is  shown  in  Figure  4.15. 
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Figure  4.15.  F-16  control  usage  with  noise,  design  point  15 

The  same  control  deflection  and  rate  problems  noted  in  Section  3.4  (for  the  same 
Wj)  are  also  evident  in  Figure  4.15.  The  current  H2/L,  methodology  could,  however,  be 
extended  to  include  Lj  control  deflection  and  rate  constraints.  This  should  be  a  topic  of 
further  research. 

4.8.  Chapter  Summary 

This  chapter  presented  a  newly  developed  numerical  approach  to  Hj/Lj 
optimization  for  continuous  systems.  This  method  was  tested  on  a  realistic  example,  and 
the  results  show  the  benefit  of  mixing  H2  and  Lj  design  problems.  The  controllers 
obtained  in  this  mixed  design  approach,  however,  must  still  be  discretized  before  they  can 
be  used  on  an  actual  aircraft.  Unfortunately,  the  performance  of  the  aircraft  with  the 
discretized  controller  can  be  significantly  different  than  the  results  obtained  with  the 
continuous  controller.  An  alternative,  and  potentially  better  design  approach  is  to  use 
mixed  H/fj  optimization  to  directly  obtain  a  discrete  controller.  This  method  is  presented 

in  Chapter  5. 
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V.  Mixed  Optimization 


5.1.  Chapter  Overview 

This  chapter  develops  a  fixed-order,  mixed  optimization  method  to  design 
controllers  for  discrete  systems.  The  chapter  is  organized  in  the  same  manner  as  Chapter 
4  to  highlight  the  similarities  and  differences  between  mixed  H2/l!i  optimization  for  discrete 
systems  and  mixed  optimization  for  continuous  systems.  Section  5.2  covers  the 
mixed  H2A1  optimization  problem.  Section  5.3  discusses  analytical  gradients  for  the  H2 
portion  of  the  problem.  Methods  of  numerically  approximating  the  norm  and  its 
gradient  are  presented  in  Sections  5.4  and  5.5,  respectively.  Computer  implementation 
issues  associated  with  optimization  are  discussed  in  Section  5.6.  In  Section  5.7,  the 
Hj/?!  optimization  method  is  demonstrated  on  the  discrete  version  of  the  AFTI  F-16 
longitudinal  model,  and  the  results  are  compared  to  the  results  obtained  in  Chapter 
4. 

5.2.  The  Mixed  Problem 

Much  of  the  mixed  optimization  problem  development  presented  in  Section 
4.2  also  applies  to  the  mixed  optimization  problem.  The  purpose  of  this  section  is  to 
simply  point  out  the  differences  between  the  two  methods. 

The  system  in  Figure  4.1,  where  the  inputs  and  outputs  are  now  sequences, 
represents  the  standard  mixed  problem.  The  goal  of  optimization  is  stated  in 
(4.1)-(4.3)  and  the  individual  H2  and  problems  can  be  expressed  in  terms  of  the  state- 
space  systems  given  in  (4.4)  and  (4.5). 

The  following  assumptions  are  made  on  the  state-space  elements  in  (4.4)  and  (4.5) 
using  the  definitions  given  in  (2.61)  and  (2.62): 

i)  (A2,  B,^)  is  stabilizable  and  (C^^,  A)  is  detectable 

ii)  RpR,>0 


5-1 


m)D^C,=0andB„D^„=0 
iv)  (Ap  B^j)  is  stabilizable  and  (C^j,  Aj)  is  detectable 
See  Section  2.9  for  the  reasons  for  assumptions  i)-iii).  Assumption  iv)  is  required  to 
ensure  the  fj  problem  has  a  solution. 

The  controller  state-space  is  given  by  (4.6),  and  the  closed-loop  state-space 
equations  for  are  given  by  (4.7).  Unlike  the  continuous  H2  problem,  a  finite  2-norm 
for  can  be  attained  when  is  non-zero  in  the  discrete  problem.  Thus,  a  strictly 
proper  controller  cannot  be  assumed  in  the  mixed  optimization  problem.  Given  this 
fact,  the  state-space  equations  for  the  mixed  H2/l!i  problem  are: 

X2  =  A2X2-hB„w 
z=C,X2-1-D,„w 

Xi  =  AjXi-l-B,r 

m=C„x,-KD^r,  (5.1) 

where 


*2  = 


(5.2) 


*1  = 


(5.3) 


A2  = 


A2  Bu2^k^y2  ^u2^k 

®k^y2  Ajj 


(5.4) 


Ai  = 


Ai  +  Bul^k^yl  ^ul^k 
®k^yl  Aj; 


(5.5) 


B,„  = 


®u2^k^yw 

^k^yw 


(5.6) 


B  = 


(5.7) 
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c.=[q  +  d^d,c,,  d^c,] 

(5.8) 

“  [^m  ^mu^k^yl  ^mu^k] 

(5.9) 

(5.10) 

Dmr  =[D„,uDkDy^-t-D„,.] 

(5.11) 

The  definitions  in  (4.18)  are  used  to  discuss  solutions  to  the  mixed  problem,  which 
must  satisfy  the  Kuhn-Tucker  necessary  conditions  given  in  Section  4.2.  An  illustration 
of  a  possible  solution  curve  is  shown  in  Figure  4.2. 

The  H2A1  problem,  like  the  H2/L,  problem,  is  solved  numerically  using  a  SQP 
algorithm.  The  objective  (f)  and  the  constraints  (g)  are  given  by 

f(c)  =  €,||T„|" 
gi(c)  =  ?.(|T.„|,-v) 

The  only  difference  between  the  equations  in  (5.12)  and  their  continuous  equivalents  in 
(4.19)  is  in  the  stability  constraint.  Recall  that  a  discrete  system  in  the  z-domain  is  stable  if 
and  only  if  its  closed-loop  eigenvalues  lie  inside  the  open  unit  disk. 

Modal  form  is  also  assumed  for  the  controller,  K,  in  the  optimization 
problem,  in  order  to  minimize  the  number  of  design  variables  in  c.  A  block- Jordan  form 
or  fully  populated  state-space  could  be  used  to  allow  for  repeated  eigenvalues  in  the 
controller. 

In  order  to  use  the  SQP  algorithm  to  numerically  solve  the  problem,  gradients 
must  be  found  for  the  objective  function  and  each  of  the  constraints.  Analytical  gradients 
for  the  stability  constraint  can  be  easily  derived  from  (4.35).  Gradients  for  one-  and  two- 
norms  are  presented  in  Sections  5.5  and  5.3,  respectively. 


5-3 


5.5.  Gradients  of  the  Two-Norm 


Recall  from  Section  2.9  that  the  two-norm  squared  of  is  given  by 

l|T„E  =  qD„Dl.+C.QCl],  (5.13) 

where  Q  is  the  positive  semidefinite  solution  to  the  Lyapunov  equation 

A,QAI-hBX=Q-  (5.14) 

If  the  constraint  on  Q  in  (5.14)  is  written  in  an  equivalent  form, 

tr[(A,QA^  +  B„B^  -Q)X]  =  0  for  all  X,  (5.15) 

then 

|T„||"  =  tr[Dl,D„  +C.QC;]+tr[(A,QAl  +B,B: -Q)X]  forall X.  (5.16) 

To  simplify  the  notation,  let  J  equal  the  two-norm  squared  of  Notice  from  (5.16) 
that  J  is  a  function  of  the  design  variables  and  the  matrices  X  and  Q.  However, 

|i=A,QAl+BX-Q  =  0.  (5.17) 

oX 

Further,  it  can  be  shown  that 


9J 

aQ 


AlXA,  +  CX-X, 


(5.18) 


which  is  a  Lyapunov  equation.  Since  X  is  arbitrary,  let  X  be  the  positive  semidefinite 
solution  to  the  Lyapunoy  equation 

A^XAj-i-C^C^-X-O.  (5.19) 

With  this  choice  of  X,  (5.18)  is  equal  to  zero.  This  choice  of  X  will  always  exist  since  A2 
must  be  stable.  The  only  remaining  derivatives  to  be  calculated  are  the  partials  of  J  with 
respect  to  the  design  variables.  If  Q  and  X  are  given  by  (4.26)  and  (4.27),  respectively, 
then  these  derivatives  are 
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(5.22) 


(5.23) 


These  equations  are  much  more  complicated  than  the  derivatives  for  the  continuous  time 
case,  because  Dj^  0. 

The  complete  method  for  finding  the  partial  derivatives  of  J  with  respect  to  the 
vectorized  compensator,  c,  is  as  follows: 

i)  Solve  the  Lyapunov  equations  in  (5.14)  and  (5.19)  for  Q  and  X,  respectively 

ii)  Compute  the  partial  derivatives  in  (5.20)-(5.23) 

iii)  Rearrange  (5.20)-(5.23)  to  form  the  gradient  with  respect  to  the  vector  of 
design  variables,  c,  as  a  vector 
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The  next  section  discusses  a  new  numerical  method  of  approximating  the  norm  of  T^. 
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5.4.  Calculating  the  fj  norm 


The  norm  of  a  discrete  SISO  system  is  given  by 


mr  I  m  i  r 
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(5.25) 


k=0 


If  Aj  is  stable,  then  A^  approaches  zero  as  k  approaches  infinity,  and  (5.25)  can  be 
approximated  by  truncating  the  infinite  series  at  some  index,  N,  which  can  be  computed 
from  the  eigenvalues  of  Aj.  The  truncated  series  can  then  be  evaluated  on  a  computer 
with  a  simple  loop,  where  k  is  the  loop  index.  With  each  pass  through  the  loop,  the 
absolute  value  of  C^A^B^  is  calculated  and  added  to  the  cumulative  sum  of  the  previous 
passes.  The  final  sum  is  added  to  the  absolute  value  of  to  produce  the  one-norm  of 
This  approach  is  used  in  the  MATLAB'''^  one-norm  calculation  routine  created  by 
Jacques  for  fixed-order,  mixed-norm  optimization  ([JR94]).  The  method  is  extremely 
reliable,  but  can  be  slow  if  N  is  large. 

The  one-norm  calculation  can  be  computed  considerably  faster,  however,  by  taking 
advantage  of  the  built-in  MATLAB'^'^  function,  Ititr.m.  This  function  rapidly  calculates  the 
vector,  (AjB^)^,  at  each  k.  The  resulting  Nxn^,  matrix  can  be  multiplied  by  to 
produce  an  Nxl  vector,  where  the  k*  element  equals  C^A^B^ .  The  absolute  sum  of 
this  vector  can  easily  be  calculated  and  added  to  the  absolute  value  of  to  produce  the 
one-norm  of  T^. 

mr 


The  norm  for  MIMO  systems  is  calculated  from  the  maximum  row  sum  of  SISO 
transfer  functions.  To  perform  this  calculation  using  this  author’s  method,  a  loop  is 
wrapped  around  the  Ititr.m  function  to  cycle  through  each  column  of  B^. 

Agorithms  based  on  both  methods  were  tested  on  a  SlUST”  SPARC  20  computer 
to  compare  efficiencies.  The  benchmark  problem  was  a  SISO  12*  order  closed-loop 
system  with  2  exogenous  inputs,  2  controlled  outputs,  and  20  design  variables  in  the 
compensator  (5*  order  compensator).  N  was  set  to  5,000.  The  first  method  took  3.5 
seconds  to  evaluate  the  one-norm  and  the  second  method  took  0.95  seconds  (all 
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performance  times  refer  to  cpu  time).  The  second  method  also  calculated  the  sign  of 
C^AjBf  for  each  column  of  and  index  k,  which  can  be  used  to  dramatically  improve 
the  speed  of  the  gradient  calculations.  In  some  parts  of  the  Yljii  algorithm  this 
computation  can  be  removed,  and  the  resulting  norm  calculation  routine  is  even  faster. 
For  example,  in  Section  5.7,  an  N  of  2  million  is  required  to  produce  an  accurate  estimate 
of  the  one-norm  for  the  AFTI  F-16  problem  at  H2  optimal.  This  computation  was 
completed  in  under  4  minutes  on  a  SUbT^  SPARC  20  computer. 

5.5.  Gradients  of  the  One-Norm 

The  analytical  gradients  of  the  one-norm  with  respect  to  the  design  variables  and 
an  improved  method  of  calculating  these  gradients  are  discussed  in  this  section.  The 
partial  derivatives  of  the  one-norm  of  T^  with  respect  to  the  closed- loop  state- space  can 
be  expressed  as 
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(5.29) 

From  these  expressions,  the  gradients  with  respect  to  the  compensator  state-space 
matrices  can  be  expressed  as 
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The  partial  of  the  one-norai  of  with  respect  to  c  can  be  formed  from  the  columns  of 
the  matrices  in  (5.30)-(5.33)  as  described  in  Section  5.3. 

Like  the  one-norm  calculation,  the  infinite  summations  in  (5.26)-(5.29)  are 
truncated  at  N,  and  can  be  calculated  with  a  series  of  loops.  This  is  the  method  used  in 
the  one-norm  gradient  routines  created  by  Jacques  for  fixed-order,  mixed-norm 
optimization  ([JR94]).  However,  large  N  values  in  this  gradient  routine  are  more  costly 
in  terms  of  computer  time  than  in  Jacques’  one-norm  calculation  routine. 

Computation  time  can  be  minimized  by  using  some  of  the  techniques  presented  in 
Section  4.5  to  compute  the  norm  for  continuous  systems.  If  Aj  is  a  non-defective 
matrix,  then  partial  of  A^  with  respect  to  any  element  of  A^  is  given  by 


^  =  J_[ra’'R-i]  =  —  A'^R-’  -t-  r|^R-'  -h  RA"  ^ .  (5.34) 
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Expressions  for  most  of  the  partial  derivatives  on  the  right  side  of  (5.34)  were  given  in 
Section  4.5.  The  only  new  derivative  to  be  calculated  is  the  partial  of  A''  with  respect  to 
ajj,  which  is  given  by 


A  3A"  dX, 


(5.35) 
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Thus,  the  partial  derivative  of  the  one-norm  with  respect  to  can  be  found  element-by¬ 
element  from 


Xsgn(C„A;B,)C„ 


k=0 
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da.:: 


(5.36) 


Recall  that  the  sign  of  C^AjB^  for  all  k  is  known  from  the  new  method  of  doing 
the  one-norm  calculation.  Therefore,  the  only  information  required  to  evaluate  (5.36)  is 
the  nj  X  Uj  diagonal  matrix.  A*" ,  at  each  index  k.  However,  since  only  the  main  diagonal 
of  this  matrix  changes,  a  Ixnj  vector  of  the  diagonal  elements  at  each  k  provides 
sufficient  information  to  reconstruct  A*" .  These  vectors  can  be  computed  rapidly  using 
Ititr.m.  The  resulting  N  x  nj  matrix  only  needs  to  be  calculated  once  in  the  gradient 
routine  regardless  of  the  number  of  inputs  and  outputs  in  the  system.  Further,  once  this 
matrix  is  calculated,  the  remaining  one-norm  derivatives  in  (5.27)-(5.33)  can  be  quickly 
evaluated. 

One-norm  gradients  for  MIMO  systems  are  calculated  for  each  output  by  summing 
the  gradients  of  the  individual  SISO  transfer  functions  between  the  output  and  the 
separate  inputs.  As  mention  in  Section  4.5,  this  approach  ensures  that  continuous  gradient 
information  is  available  regardless  of  where  the  maximum  row  sum  occurs. 

Both  methods  of  calculating  the  one-norm  gradients  were  tested  on  the  benchmark 
problem  described  in  the  last  section.  The  first  method  took  28.2  hours  to  complete  one 
gradient  calculation  on  a  SUbT'^  SPARC  20  computer.  The  second  method  took  5.03 
seconds.  It  should  be  noted  that  the  first  method  works  relatively  fast  for  systems  where 
N  is  small.  However,  larger  N  values  are  required  for  systems  with  both  fast  and  slow 
dynamics.  For  example,  the  AFTI F-16  design  problem  presented  in  Section  5.7  required 
N  values  of  at  least  5,000. 

The  new  approach  to  the  one-norm  gradient  calculation  requires  that  the 
eigenvalues  of  Aj  be  distinct.  If  this  condition  is  violated,  the  gradient  routine  switches  to 
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finite  difference  calculations  based  on  the  new  approach  of  computing  the  norm.  While 
less  accurate  than  the  old  method  of  computing  the  gradients  of  the  one-norm,  the  finite 
difference  calculations  are  still  considerably  faster. 

5.6.  Computer  Implementation 

The  author’s  norm  and  gradient  routines  were  added  to  the  fixed-order,  mixed- 
norm  optimization  algorithm  created  by  Jacques  ([JR94]).  This  algorithm  is  capable  of 
mnning  several  different  types  of  discrete  optimization  problems  with  different  SQP 
routines.  Section  5.7  demonstrates  the  H2/l!i  portion  of  the  algorithm  using  the 
MATLAB’’’^  SQP  routine,  constr.m. 

5 . 7.  Design  Example 

The  optimization  method  is  demonstrated  on  a  discrete  version  of  the  H2/L] 
design  problem  presented  in  Section  4.7.  The  state-space  matrices  for  the  discrete  F-16 
plant  are  given  in  Appendix  A.  The  discrete  state-space  matrices  for  the  H2  and  sub¬ 
problems  are  shown  in  Appendix  C. 

5.7.1.  Results 

A  plot  of  the  HJk  solution  curve  for  a  5th  order  compensator  is  shown  in  Figure 

5.1.  Values  for  several  of  the  solution  points  are  shown  in  Table  5.1.  Point  number  17 
represents  the  discrete  H2  optimal  solution,  which  is  not  depicted  in  Figure  5.1.  The 
points  in  Table  5.1  are  referred  to  as  the  discrete  design  points  for  the  remainder  of  this 
section. 
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Plots  of  the  sensitivity  and  complementary  sensitivity  for  the  discrete  design  points 


are  shown  in  Figures  5.2  and  5.3. 


Figure  5.3.  F-16  complementary  sensitivity  from  discrete  design 


The  distinctively  different  curve  in  each  plot  corresponds  to  the  H2  solution. 
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The  vector  gain  and  phase  margins  for  the  discrete  design  points  are  shown  in 


Table  5.2. 


point  # 

lower  VGM  (dB) 

upper  VGM  (dB) 

VPM  (deg) 

1 

-11.2 

4.73 

42.4 

4 

-12.6 

6.70 

45.0 

7 

-9.80 

9.73 

39.5 

11 

-8.01 

10.8 

41.7 

14 

-7.61 

9.97 

39.9 

17 

-10.5 

11.4 

42.9 

Table  5.2.  F-16  vector  gain  and  phase  margins  from  discrete  design 


Notice  that  the  stability  margins  do  not  consistently  improve  as  v  decreases,  but  do  remain 
acceptable  for  all  the  discrete  design  points. 

The  open-loop  GK  plot  is  shown  in  Figure  5.4 
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Notice  that  the  mixed  designs  miss  the  performance  and  disturbance  rejection 
barrier.  The  sensor  noise  and  unmodeled  dynamics  barrier  is  not  depicted  since  digital 
aliasing  would  cause  it  to  be  violated.  An  anti-aliasing  filter  (basically  a  low-pass  filter) 
can  be  added  to  the  system  to  produce  any  amount  of  high  frequency  roll-off. 

Figure  5.5  shows  the  F-16  responses,  without  noise,  to  a  commanded  Ig  step  input 
for  the  different  designs.  Notice  that  the  mixed  designs  track  with  no  steady-state  error. 
Step  responses  with  noise  for  discrete  design  points  1,7,  14  and  17  are  shown  in  Figures 
5. 6-5.9.  These  plots  demonstrate  that  tracking  performance  can  be  greatly  improved 
using  mixed  Hj/iSi  designs  with  very  little  increase  of  noise  in  the  system  response.  High 
frequency  noise  is  not  as  prevalent  in  these  plots  as  it  was  in  the  continuous  simulations 
due  to  the  sampling  of  the  discrete  controller. 


Figure  5.5.  F-16  step  responses  without  noise,  from  discrete  design 
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Figure  5.6.  F-16  step  response  with  noise,  discrete  design  point  1 


Figure  5.7.  F-16  step  response  with  noise,  discrete  design  point  7 


5-15 


Figure  5.8.  F-16  step  response  with  noise,  discrete  design  point  14 
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Discrete  design  point  7  is  chosen  as  the  best  mixed  H2A1  design  based  on  stability 
margins,  loop  shape  and  step  responses.  The  control  usage  for  this  particular  discrete 
design  point  is  shown  in  Figure  5.10.  Discrete  design  point  7  is  compared  to  the  best 
Hj/Lj  results  from  Chapter  4  in  the  next  section. 


Figure  5.10.  F-16  control  usage  with  noise,  discrete  design  point  7 

5.7.2.  Comparing  HJL.^  to 

It  is  impossible  to  pick  an  equivalent  point  from  the  mixed  Hj/Lj  solutions  and  the 
mixed  solutions  for  direct  comparison.  However,  some  general  conclusions  can  be 

drawn  about  the  two  methods.  To  illustrate  the  points  in  this  section,  the  AFTI  F-16  is 
tested  with  the  best  controller  discretized  at  30Hz  using  a  ZOH.  The  step  response 
and  control  usage  are  shown  in  Figures  5.1 1  and  5.12,  respectively. 
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Notice  that  the  discretized  controller  does  not  produce  the  same  level  of 
performance  as  it  did  prior  to  discretization  (see  Figures  4.13  and  4.15).  The  system  with 
the  discretized  controller  has  larger  responses  to  noise. 

The  responses  in  Figures  5.11  and  5.12  are  similar  to  the  responses  for  the  best 
H/j  controller  (Figures  5.7  and  5.10)  in  terms  of  peak  values,  peak  times  and  settling 

times.  This  suggests  that  comparable  designs  could  be  done  with  either  method. 
However,  the  Hj/Lj  controller  is  6th  order  and  the  controller  is  5th  order.  If 
compensator  orders  are  chosen  to  be  the  size  of  the  H2  plant  in  either  design,  then  IIJL^ 
optimization  will  always  produce  a  higher  order  compensator.  This  is  due  to  the 
additional  dynamics  that  must  added  to  the  continuous  design  to  represent  the  time  delay 
that  will  occur  when  the  controller  is  discretized  (see  Appendix  A). 

The  Fade  approximations  used  to  represent  the  time  delay  in  the  design  can 
cause  additional  problems.  First,  they  make  the  continuous  system  nonminimum  phase. 
This  can  be  seen  in  the  Hj/Lj  step  responses  in  Figure  4.10,  where  all  initially  respond  in 
the  opposite  direction  of  the  commanded  normal  acceleration.  Since  the  actual  aircraft 
model  is  minimum  phase  without  the  Fade  approximations,  these  plots  can  be  misleading. 
Second,  large  order  Fade  approximations  are  required  to  accurately  represent  a  time 
delay,  which  can  lead  to  higher  order  controllers  in  an  H^/Lj  design.  If  a  lower  order  Fade 
approximation  is  used,  care  must  taken  to  discretize  the  continuous  controller  at  a  slighdy 
faster  period  than  the  delay  time  represented  by  the  approximation.  Discretizing  the 
controller  at  a  period  equal  to  the  delay  time  represented  by  the  approximation  can  lead  to 
an  unstable  closed-loop  system.  In  light  of  the  issues  covered  above,  it  can  be  concluded 
that  digital  controller  designs  are  better  accomplished  with  H2/{i  optimization. 

5.8.  Chapter  Summary 

This  chapter  presented  a  numerical  approach  to  optimization  for  discrete 
systems.  New  and  efficient  methods  of  approximating  the  one-norm  and  its  gradients 
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were  also  presented.  The  H2/l!i  algorithm  was  tested  on  a  realistic  problem  which  showed 
the  benefits  of  this  mixed-norm  design  approach.  Comparisons  of  systems  with  discrete 
BySj  controllers  and  discretized  controllers  showed  that  adequate  digital  controllers 
are  better  obtained  with  HJl^  optimization.  Chapter  6  summarizes  the  results  in  this  thesis 
and  presents  recommendations  for  further  research. 
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VI.  Conclusions  and  Recommendations 


6.1.  Summary 

This  thesis  explored  the  use  of  optimization  for  flight  controller  designs. 
Methods  of  using  optimization  were  developed  to  handle  “hard”  magnitude  constraints, 
such  as  control  deflection  limitations,  control  rate  limitations  and  overshoot  restrictions 
which  are  difficult  to  incorporate  into  other  optimization  methods.  These  constraints  were 
added  to  the  standard  ffee-order  ly  optimization  problem  and  demonstrated  on  a  realistic 
flight  control  design  with  mixed  success.  The  ly  optimization  algorithm  worked  well  for 
small  systems  with  one  or  two  constraints,  but  had  difficulty  handling  very  large  systems 
of  equations  with  multiple  constraints.  The  constrained  ly  optimization  problem  also  had 
the  tendency  to  produce  compensators  with  unreasonably  large  orders.  To  counter  this 
problem,  two  model  matching  techniques  were  presented  based  on  one  block  systems. 
These  one  block  problems  were  solved  quickly  and  accurately  with  the  current  ty 
optimization  software,  and  acceptable  tracking  results  were  produced  for  SISO  and 
MIMO  systems  with  controller  orders  less  than  or  equal  to  the  weighted  plant. 

Hj  optimization  was  mixed  with  ty  optimization  to  produce  systems  with  good 
noise  performance  and  tracking  characteristics.  This  method  was  done  first  in  continuous 
time  with  fixed-order  H^/Ly  optimization.  New  numerical  approaches  of  calculating  the  Lj 
norm  and  its  gradient  were  developed  and  incorporated  into  the  iiy/Ly  algorithm  which 
avoid  many  of  the  complications  associated  with  previous  methods  dependent  on  the 
forward  Euler  transformation.  The  K/Ly  algorithm  was  tested  on  a  SISO  system  and  the 
results  showed  the  benefits  of  a  mixed  design  over  a  pure  H2  or  Ly  design.  This  method 
had  the  added  benefit  that  the  compensator  order  is  completely  chosen  by  the  designer. 

Fixed-order  YLJly  optimization  was  then  used  to  design  discrete  systems.  New 
numerical  approaches  of  calculating  the  ly  norm  and  its  gradient  were  presented  which 
offered  dramatic  improvements  in  terms  of  computational  efficiency,  and  essentially 
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extended  the  utility  of  the  existing  H^/IEi  algorithm  to  a  wider  class  of  problems.  The 
algorithm  was  tested  on  a  realistic  SISO  system  that  could  not  have  been  solved  prior  to 
the  improvement  in  the  algorithm.  The  results  showed  that  a  good  blend  of  tracking  and 
noise  performance  was  easy  to  obtain  with  mixed  optimization.  Like  the  Hj/L, 
method,  the  size  of  the  controller  in  optimization  is  chosen  by  the  designer. 

Two  different  methods  of  digital  controller  designs  were  examined  in  the  context 
of  Hj/Lj  optimization  and  optimization.  The  first  method  involved  completing  a 
satisfactory  Hj/Lj  design  and  discretizing  the  resulting  controller.  In  the  second  method, 
the  discrete  controller  was  obtained  directly  with  H2/i!i  optimization.  Comparisons  of  the 
two  approaches  showed  that  comparable  designs  could  be  done  with  either  method,  but 
the  Hj/?!  approach  would  generally  produce  a  smaller  order  controller  due  to  the  additional 
dynamics  required  in  the  problem  to  approximate  a  time  delay.  Since  the 
algorithm  is  also  quicker  and  more  accurate  than  the  Hj/Lj  algorithm,  it  can  be  concluded 
that  fixed-order  digital  controller  design  is  best  done  with  HjAi  optimization. 

Of  all  the  methods  presented  in  this  thesis,  fixed-order  optimization  offers  the 
most  promise  for  successful  designs  of  small  order  digital  flight  controllers  that  produce 
systems  with  good  noise  and  tracking  performance.  However,  research  still  needs  to  be 
done  on  free-order  optimization  to  determine  the  amount  of  performance  lost  in 
fixed-order  optimization  by  setting  the  size  of  the  controller  a  priori.  Additional 
research  topics  are  discussed  in  the  next  section. 

6.2 .  Recommendations  for  Future  Research 

While  this  thesis  has  contributed  to  a  better  understanding  of  how  and  mixed 
Hj/^i  optimization  can  be  applied,  more  research  needs  to  be  done  in  these  two  areas. 
First,  the  optimization  algorithm  needs  to  be  improved  to  allow  it  to  handle  large 
systems  with  multiple  constraints.  This  could  probably  be  accomplished  with  an  improved 
linear  programming  algorithm  and  appropriate  scaling  of  the  design  variables.  Additional 
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improvements  are  required  to  increase  the  speed  of  the  algorithm,  and  to  simplify  the 
way  constraints  are  augmented  to  a  standard  design  problem.  Further  research  in 
frequency  constraints  (not  discussed  in  this  thesis)  should  also  be  conducted. 

Second,  many  of  the  “hard”  constraints  incorporated  into  optimization  should  be 
incorporated  into  optimization.  This  should  not  be  difficult  to  do  given  the  present 
structure  of  the  U^/li  algorithm.  The  added  capabilities,  however,  will  greatly  improve  the 
utility  of  the  portion  of  the  U^/l^  algorithm. 

Third,  the  H^/l^  algorithm  requires  several  modifications  before  it  can  be  reliably 
used  on  wide  range  of  problems.  Research  needs  to  be  done  to  allow  the  analytical 
gradients  to  handle  systems  with  repeated  eigenvalues.  The  current  version  of  the 
software  relies  on  finite  difference  calculations  in  this  case.  Repeated  eigenvalues  are  also 
not  permitted  in  the  resulting  controller.  This  problem  could  be  fixed  by  allowing  a 
Jordan  form  for  the  controller  instead  of  restricting  it  to  modal  form.  Further,  scaling 
problems  in  first  and  second  order  derivatives  in  the  MATLAB’’’^  SQP  routine,  constr.m, 
restrict  the  ability  of  the  HJli  algorithm  to  quickly  generate  a  complete  design  curve  (such 
as  the  one  depicted  in  Figure  5.1).  While  the  IMSL''’^  SQP  routine  alleviates  some  of 
these  problems,  it  is  commercial  software  and  thus  can  not  be  modified  to  suit  specific 
requirements. 

Finally,  further  research  is  needed  on  mixed  optimization  for  discrete 

systems.  The  portion  of  the  problem  could  be  used  to  improve  stability  margins  with 
sensitivity  minimization.  This  capability  was  not  possible  with  sensitivity  minimization. 
The  portion  of  the  problem  could  then  be  used  to  handle  the  “hard”  magnitude 
constraints  of  the  system.  Work  is  currently  underway  in  this  area. 
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Appendix  A.  SISO  AFTI F-16  Model 


The  AFTI  F-16  model  for  normal  acceleration  command  following  is  shown  in 
Figure  A.l.  The  continuous  system  consists  of  a  4  state  model  of  the  aircraft’s 
longitudinal  dynamics  (Wp),  a  first  order  actuator  (WJ  and  a  first  order  Fade 
approximation  of  a  0.05  second  time  delay  (W^).  The  time  delay  model  is  not  required  in 
the  discrete  version  of  the  system. 


Wi 


Figure  A.l.  F-16  model  block  diagram 


The  four  states  in  the  longitudinal  model  are  forward  speed  (u  in  ft/sec),  angle  of 
attack  (a  in  radians) ,  pitch  angle  (0  in  radians),  and  pitch  rate  (q  in  radians/sec).  The 

input  to  Wp  is  the  stabilator  deflection  (6^  in  radians)  and  the  output  is  the  normal 
acceleration  (n^  in  g’s).  Wp  is  given  by 


u 

■-1.485e-2 

3.738e  +  l 

-3.220e  +  l 

-1.794e-H' 

u 

■  2.140e-3‘ 

d 

-8.000e-5 

-1.491e  +  0 

-1.300e-3 

9.960e-l 

a 

1 

-1.880e-l 

0 

O.OOOe  +  O 

O.OOOe  +  0 

O.OOOe  +  0 

1.000e-F0 

0 

T 

O.OOOe  +  0 

q. 

-3.600e-4 

9.753e  +  0 

2.900e-4 

-9.600e-l 

.q. 

_-1.904e-i-l_ 

A-1 


n,  =[l.500e-3  35264e  +  l  2.720e-2  -3.340e-l] 


+  [-4.366e  +  0]6,. 


The  input  to  is  the  commanded  stabilator  deflection  ( 5^^  in  radians)  and  the 
output  is  the  stabilator  deflection.  W,  is  given  by 


X,  =  [-2.000e  +  l]x,  +  [2.000e  + 1]5,^ 

5,  =  [LOOOe  +  0]x,  +  [O.OOOe  +  0]5,^ . 

The  input  to  is  normal  acceleration  and  the  output  is  the  delayed  normal 
acceleration  ( n^^  in  g’s  ).  is  given  by 


Xj  =  [-4.000e  +  l]xj  +  [l.000e  +  0]n^ 
n^^  =  [S.OOOe  +  l]x,  +  [— l.OOOe  +  0]n^ . 

The  continuous  aircraft  plant,  G^,  equals  W^WpWj.  G^.  is  given  by 


G 


c 


Ac 

Be 

.Cc 

Do 

where 


A. 


-1.485e-2 
-S.OOOe -5 
O.OOOe +  0 
-3.600e-4 
O.OOOe +  0 
1.500e-3 


3.738e  +  l 
-1.491e  +  0 
O.OOOe  +  0 
9.753e  +  0 
O.OOOe  +  0 
3.526e+l 


-3.220e+l 
-1.300e-3 
O.OOOe +  0 
2.900e-4 
O.OOOe +  0 
2.720e-2 


-1.794e  +  l 
9.960e-l 
LOOOe +  0 
-9.600e-l 
O.OOOe +  0 
-3.340e-l 


2.140e-3 
-1.880e-l 
O.OOOe  +  0 
-1.904e  +  l 
-2.000e  +  l 
-4.366e  +  0 


O.OOOe +  0 
O.OOOe +  0 
O.OOOe +  0 
O.OOOe +  0 
O.OOOe +  0 
-4.000e  +  0 
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B 


C 


O.OOOe  +  0 
O.OOOe  +  0 
O.OOOe  +  0 
O.OOOe  +  0 
2.000e  +  l 
O.OOOe  +  0 


C,  =  [  -liOOe  -  3  -3^26e  + 1  -2.720e  -  2  3.340e  - 1  4.366e  +  0  S.OOOe  + 1  ] 

D,  =[0.000e  +  0]. 


The  discrete  aircraft  plant,  G^,  equals  W^Wp  discretized  at  30Hz  using  a  ZOH.  Gj 
is  given  by 


where 


9.995e-l 
-2.800e-6 
-2.028e-7 
-1.225e-5 
O.OOOe  +  0 


1.121e  +  0 
9i67e-l 
5.278e-3 
3.127e-l 
O.OOOe  +  0 


-1.073e  +  0 
-4.072e-5 
l.OOOe  +  0 
9.198e-6 
O.OOOe  +  0 


-5.870e-l 
3.193e-2 
3.286e-2 
9.737e-l 
O.OOOe  +  0 


1.489e-l 

-1.276e-2 

-8.494e-3 

-4.568e-l 

5.134e-l_ 


C, 


B 


d 


3.494e-2 

-3.623e-3 

-1.992e-3 

-1.700e-l 

4.866e-l 


[l.500e-3  3i26e+l  2.720e-2  -3.340e-l  ^.366e  +  0] 


D,  =[0.000e  +  0]. 


The  truth  model  used  in  simulations  includes  a  Von  Karman  wind  model  (WJ  for 
the  wind  noise,  and  a  high-pass  filter  (W  J  for  the  measurement  noise.  is  given  by 
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X  „  =  [-6.700e  +  0]x  ^  +  [l.870e  -  3]w  j 
^i=[l.000e  +  0]x„. 

where  Wj  is  a  unit  strength  white  Gaussian  noise,  and  is  the  wind  noise.  is  given  by 

x^  =  [-l.OOOe  +  l]x^  +  [4.000e  -  4]w2 
^2  =  [LOOOe  +  0]x„  +  [4.000e  -  4]w2 . 

where  W2  is  a  unit  strength  white  Gaussian  noise,  and  ^2  is  the  measurement  noise.  The 
wind  noise  enters  the  aircraft  plant  as  an  angle  of  attack  perturbation.  F,  in  Figure  A.l, 
equals  the  second  column  of  the  A  matrix  in  the  longitudinal  model. 
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Appendix  B.  MIMO  Missile  Model 


The  3  input,  3  output  missile  model  (taken  from  [Bro91]  and  [Luk93])  for  normal 
and  lateral  acceleration  command  following  is  shown  in  Figure  B.l.  The  system  consists 
of  a  5  state  model  of  the  missile’s  lateral  and  longitudinal  dynamics  (Wp),  and  a  3  state 
model  of  the  actuator  dynamics  (W^). 


Figure  B.l.  Missile  model  block  diagram 


The  5  states  in  Wp  are  angle  of  attack  (a  in  radians) ,  sideslip  angle  (P  in  radians), 
roll  rate  (p  radians/sec),  pitch  rate  (q  radians/sec),  and  yaw  rate  (r  radians/sec).  The  inputs 
to  Wp  are  the  control  fin  deflections  (5p,  8^,  8,  in  radians)  and  the  outputs  are  the  roll  rate, 
lateral  acceleration  (n^  in  ft/sec^),  and  normal  acceleration  (n^  in  ft/sec^).  Wp  is  given  by 


-9.844e-l  -9.234e-2  O.OOOe  +  0  l.OOOe  +  0 

-9.234e-2  -9.844e-l  O.OOOe  +  0  O.OOOe  +  0 

-2.674e  +  2  2.674e  +  2  -1.595e  +  0  O.OOOe  +  0 

-1.946e  +  2  -6.967e-l  O.OOOe  +  0  -1550e  +  0 

6.967e-l  1.946e  +  2  O.OOOe  +  0  O.OOOe  +  0 

■-2.186e-3  3.26  le-1  -1.135e-3' 

-6.189e-4  1.566e-3  -3.188e-l  fdp' 

+  1.140e  +  4  -6.550e  +  l  -1.070e  +  l  8, 

-2.736e  +  0  3.21  le  + 2  1.451e-l  [8,_ 

-3.045e  - 1  1.4  lOe  +  0  3.149e  +  2 


O.OOOe  +  0  a 
-l.OOOe  +  0  p 
O.OOOe  +  0  p 
O.OOOe +  0  q 
-li50e  +  0  r 
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O.OOOe  +  0 

aoooe+0 

l.OOOe  +  0 

O.OOOe  +  0 

O.OOOe  +  0 

— 2.062e  +  2 

-2.198e  +  3 

O.OOOe  +  0 

O.OOOe  +  0 

O.OOOe  +  0 

-2.198e  +  3 

-2.062e  +  2 

O.OOOe  +  0 

O.OOOe  +  0 

O.OOOe  +  0 

■  O.OOOe  +  0 

O.OOOe  +  0 

0.000e  +  0‘ 

O.OOOe  +  0 

O.OOOe  +  0 

O.OOOe  +  0 

■Sp' 

O.OOOe  +  0 

O.OOOe  +  0 

O.OOOe  +  0 

Sq 

-1.368e  +  0 

3.496e  +  0 

-7.119e  +  2 

S.. 

-4.876e  +  0 

7.282e  +  2 

-2J35e  +  0 

The  inputs  to  are  the  commanded  fin  deflections  (5p^,  8^^,  8^^  in  radians)  and 
the  outputs  are  the  actual  fin  deflections.  is  given  by 


32 


-4.000e  +  2  O.OOOe  +  0  O.OOOe  +  0 

O.OOOe  +  0  -4.000e  +  2  O.OOOe  +  0 


*31 


32 


L 

O.OOOe  +  0 

O.OOOe  +  0 

-4.000e  +  2j[x 

"LGOOe  +  O 

O.OOOe  +  0 

0.000e  +  0l 

+ 

O.OOOe  +  0 

l.OOOe  +  0 

O.OOOe  +  0 

5. 

O.OOOe  +  0 

O.OOOe  +  0 

1.000e  +  oJ 

.5j 

■5p' 

4.0006  +  2 

O.OOOe  +  0 

O.OOOe  +  0 

X' 

5q 

= 

O.OOOe  +  0 

4.000e  +  2 

O.OOOe  +  0 

^32 

6. 

O.OOOe  +  0 

O.OOOe  +  0 

4.000e  +  2_ 

.^33. 

rO.OOOe  +  0 

O.OOOe  +  0 

O.OOOe  +  0 

+ 


O.OOOe  +  0  O.OOOe  +  0  O.OOOe  +  0 
O.OOOe  +  0  O.OOOe  +  0  O.OOOe  +  0 


qc 
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The  discrete  aircraft  plant,  G^,  equals  W^W  discretized  at  lOOHz  using  a  ZOH. 


is  given  by 
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1.832e-2 
O.OOOe  +  0 
O.OOOe  +  0 
-2.250e-2 
1.679e-3 
1.105e  +  4 
-2.634e  +  0 
-2.945e-l 
O.OOOe  +  0 
O.OOOe  +  0 
O.OOOe  +  0 
-9.428e-4 
9.806e-l 
2.632e  +  0 
-5.950e  +  0 
1.916e  +  0 


O.OOOe  +  0 
1.832e-2 
O.OOOe  +  0 
2.708e  +  0 
-1.015e-2 
-6.682e  +  l 
3.09  le  + 2 
1.369e  +  0 
O.OOOe  +  0 
O.OOOe  +  0 
O.OOOe  +  0 
O.OOOe  +  0 
O.OOOe  +  0 
9.842e-l 
O.OOOe  +  0 
O.OOOe  +  0 


O.OOOe  +  0 
O.OOOe  +  0 
1.832e-2 
U21e-3 
-2.655e  +  0 
-1.360e  +  l 
1.492e-l 
3.032e  +  2 
O.OOOe  +  0 
O.OOOe  +  0 
O.OOOe  +  0 
9.842e-3 
-4.663e-6 
-1.317e-2 
9.750e-l 
3.130e-5 


O.OOOe  +  0 
O.OOOe  +  0 
O.OOOe  +  0 
9.806e-l 
-9.428e-4 
-2.632e  +  0 
-1.916e  +  0 
5.950e-3 
O.OOOe  +  0' 
O.OOOe  +  0 
O.OOOe  +  0 
4.663e-6 
-9.842e-3 
-1.317e-2 
3.130e-5 
9.750e-l 


2.4542e-3 

O.OOOe +  0 

0.000e 

O.OOOe  +  0 

2.4542e-3 

O.OOOe 

O.OOOe  +  0 

O.OOOe +  0 

2.4542e 

-1.009e-4 

1.236e-2 

-5.034e 

4.790e-6 

-3.539e-5 

-1.21  le 

8.54  le  +  1 

-5.013e-l 

-9.025e 

-2.045e-2 

2.400e  +  0 

1.1 17e 

-2.282e-3 

1.058e-2 

2.354e 

O.OOOe  +  0 
-5.474e  +  2 
-1.9506  +  3 
O.OOOe  +  0 
-2.198e  +  3 
-2.062e  +  2 


O.OOOe  +  0 
1.398e  +  3 
2.913e  +  5 
l.OOOe  +  0 
O.OOOe  +  0 
O.OOOe  +  0 


O.OOOe  +  0 
-2.848e  +  5 
-1.014e  +  3 

O.OOOe  +  0 
O.OOOe  +  0 
O.OOOe  +  0 


O.OOOe  +  0 
— 2.062e  +  2 
-2.19796  +  3 
0.0006  +  0' 
0.0006  +  0 
0.0006  +  0 


O.OOOe  +  0  0.0006  +  0  0.0006  +  0 
0.0006  +  0  0.0006  +  0  0.0006  +  0 


Dd 


0.0006  +  0  0.0006  +  0 
0.0006  +  0  0.0006  +  0 
0.0006  +  0  0.0006  +  0 


0.0006  +  0 
0.0006  +  0 
0.0006  +  0 


Appendix  C.  Matrices  for  and  Subproblems 


The  matrices  for  the  continuous  H2  problem  are: 


■-1.485e-2 

3.738e  +  l 

-3.220e  +  l 

-1.794e  +  l 

2.140e-3 

-8.000e-5 

-1.49  le  +  0 

-1.300e-3 

9.960e-l 

-1.880e-l 

O.OOOe  +  0 

O.OOOe  +  0 

O.OOOe  +  O 

l.OOOe  +  0 

O.OOOe  +  0 

A,  = 

2 

-3.600e-4 

9.753e  +  0 

2.900e-4 

-9.600e-l 

-1.904e  +  l 

O.OOOe  +  0 

O.OOOe  +  0 

O.OOOe  +  0 

O.OOOe  +  0  - 

-2.000e  +  l 

1.500e-3 

3i26e+l 

2.720e-2 

-3.340e-l  - 

-4.366e  +  0 

8.359e-l 

O.OOOe  +  0' 

"0.0006  + O' 

-3.334e-2 

0.0006 +  0 

O.OOOe  +  0 

O.OOOe  +  0 

O.OOOe  +  0 

O.OOOe  +  0 

B  = 

Bu2  = 

2.181e-l 

O.OOOe  +  0 

O.OOOe  +  0 

O.OOOe  +  0 

O.OOOe  +  0 

2.000e  +  l 

O.OOOe  +  0 

O.OOOe+0 

O.OOOe  +  0 

■-l^OOe-3 

-3.526e  +  l 

-2.720e-2 

3.340e-l 

4.366e  +  0 

C,  = 

O.OOOe  +  0 

z 

O.OOOe  +  0 

O.OOOe  +  0 

O.OOOe  +  0 

O.OOOe  +  0 

[-1500e-3 

-3i26e+l 

-2.720e-2 

3.340e-l 

4.366e  +  0 

D, 


"0.0006  +  0  O.OOOe  +  0" 

D™  = 

"0.0006  + 0' 

O.OOOe  +  0  O.OOOe  +  0 

zu 

1.000e  +  l_ 

Dy,  =[0.000e  +  0  4.000e-3] 


D^=[0.000e  +  0] 


O.OOOe  +  0 
O.OOOe  +  0 
O.OOOe  +  0 
O.OOOe  +  0 
O.OOOe  +  0 
-4.000e+l 


8.000e  +  l 
O.OOOe  +  0 

8.000e  +  l] 
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The  matrices  for  the  Lj  problem  are 


1 

oo 

1 

to 

3.738e  +  l 

-3.220e  +  l 

-8.000e-5 

-L491e  +  0 

-L300e-3 

O.OOOe  +  0 

O.OOOe  +  0 

O.OOOe +  0 

-3.600e-4 

9.753e  +  0 

2.900e-4 

O.OOOe +  0 

O.OOOe +  0 

O.OOOe +  0 

1300e-3 

3526e  +  l 

2.720e-2 

-1300e-3 

-3326e  +  l 

-2.720e-2 

... 

-1.794e  +  l 

2.140e-3 

O.OOOe +  0 

O.OOOe +  0 

9.960e-l 

1 

00 

00 

T 

O.OOOe +  0 

O.OOOe +  0 

LOOOe +  0 

O.OOOe +  0 

O.OOOe +  0 

O.OOOe +  0 

-9.600e-l 

-L904e  +  1 

O.OOOe +  0 

O.OOOe +  0 

O.OOOe +  0 

-2.000e  +  l 

O.OOOe +  0 

O.OOOe +  0 

-3.340e-l 

-4.366e  +  0 

-4.000e+l 

O.OOOe +  0 

3.340e-l 

4.366e  4-  0 

8.000e  +  l 

-LOOOe -4 

O.OOOe  +  O' 

‘O.OOOe  +  0' 

O.OOOe +  0 

O.OOOe  +  0 

O.OOOe +  0 

O.OOOe  +  0 

O.OOOe +  0 

O.OOOe  +  0 

O.OOOe +  0 

2.000e  +  l 

O.OOOe +  0 

O.OOOe  +  0 

LOOOe +  0 

O.OOOe  +  0 

C^=[-U00e-3  -3526e  +  l  -2.720e-2 

3.340e-l  4366e  +  0  8.000e  +  l  1.000e  +  l] 
Cyi  =[-1300e-3  -3326e  +  l  -2.720e-2  ••• 

3.340e-l  4.366e  +  0  8.000e  +  l  O.OOOe  +  O] 

=  [LOOOe  +  0]  =  [O.OOOe  +  O] 

D^=  [LOOOe  +  0]  =  [O.OOOe  +  O] 

C-2 


The  state-space  matrices  for  the  discrete  H2  problem  are 


■  9.995  le-1 

1.1216  4-0 

-1.073e-h0 

-5.870e-l 

1.489e-l 

-2.800e-6 

9.567e-l 

-4.072e-5 

3.193e-2 

-1.276e-2 

-2.028e-7 

5.278e-3 

LOOOe-^O 

3.286e-2 

-8.494e-3 

-1.225e-5 

3.127e-l 

9.198e-6 

9.737e-l 

-4568e-l 

O.OOOe  +  0 

0.0006  4-0 

0.0006  4-0 

O.OOOe  +  0 

0.0006  4-0 

■  2.506e-2 

O.OOOe-i-O' 

'  3.494e-2' 

-9.686e-4 

0.0006  4-0 

-3.623e-3 

1.180e-4 

0.000e4-0 

-1.992e-3 

6.99  le- 3 

O.OOOe-i-0 

-1.699e-l 

O.OOOe  +  0 

0.000e4-0 

_  4.865e-l_ 

C 


z 


-1.500e-3 
O.OOOe  +  0 


-3.564e-M  -2.720e-2 
O.OOOe  +  0  O.OOOe-t-0 


Cy2  =  [  -liOOe  -  3  -3^64e  1  -2.720e  -  2 


3.340e-l  4.366e  +  0 
O.OOOe-t-0  O.OOOe  +  0 
3.340e-l  4.366e-H0] 


"0.0006  4-0 

0.000e4-0" 

= 

"O.OOOe-i-O" 

0.0006  4-0 

0.0006  4-0 

zu 

_1.000e-H_ 

[O.OOOe  +  0 

4.000e-3] 

[O.OOOe-HO] 

C-3 


The  matrices  for  the  problem  are 


9.995  le-1 

1.121e  +  0 

-1.073e  +  0 

-5.870e-l 

1.489e-l 

-2.800e-6 

9i67e-l 

-4.072e-5 

3.193e-2 

-1.276e-2 

-2.028e-7 

5.278e-3 

l.OOOe  +  0 

3.286e-2 

-8.494e-3 

-1.225e-5 

3.127e-l 

9.198e-6 

9.737e-l 

-4.568e-l 

O.OOOe  +  0 

O.OOOe  +  0 

O.OOOe +  0 

O.OOOe +  0 

5.134e-l 

-4.805e-5 

-1.147e  +  0 

-8.555e-4 

-8.030e-3 

1.098e-l 

'O.OOOe  +  O' 

■  3.494e-2' 

O.OOOe +  0 

-3.623e-3 

O.OOOe +  0 

= 

-1.992e-3 

O.OOOe +  0 

ul 

-1.699e-l 

O.OOOe +  0 

4.865e-l 

_3.333e-2_ 

_  3.994e-2_ 

=  [  -  1.500e  -  3  -3^64e  + 1  -2.720e  -  2 

Cyi  =  [  -1.500e  -  3  -3.564e  + 1  -2.720e  -  2 


3.340e-l  4.366e  +  0 
3.340e-l  4.366e  +  0 


D..=[l.000e  +  0]  D^„=[0.000e  +  0] 

D^=  [l.000e  +  0]  D^=  [O.OOOe  +  O] 


C-4 


O.OOOe +  0 
O.OOOe +  0 
O.OOOe +  0 
O.OOOe +  0 
O.OOOe +  0 
l.OOOe  +  0 


l.OOOe  +  1] 
O.OOOe  +  0] 
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